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Gauss  Elimination: 

Workhorse  of  Linear  Algebra 

Peter  R  Turner 

Mathematics  Department,  U  S  Naval  Academy,  Annapolis,  MD  21402 

Abstract.  This  report  brings  together  many  different  aspects  of  Gauss  elimi¬ 
nation.  The  basic  Gauss  elimination  (GE)  algorithm  is  a  fundamental  tool  of  linear 
algebra  computation  for  solving  s>'stems,  computing  determinants  and  determining 
the  rank  of  a  matrix.  All  of  these  are  discussed  in  varying  contexts.  These  include 
different  arithmetic  or  algebraic  settings  such  as  integer  arithmetic  or  polynomial 
rings  as  well  as  conventional  real  (floating-point)  arithmetic.  These  have  effects  on 
both  accuracy  and  complexity  analyses  of  the  algorithm.  These,  too,  are  covered 
here.  The  impact  of  modern  parallel  computer  architecture  on  GE  is  also  included. 
Finally,  GE  is  considered  within  the  context  of  “noisy”  matrices.  The  effect  of  the 
noise  in  matrix  entries  on  the  effective  rank  of  the  matrix  is  the  central  aspect  con¬ 
sidered  here. 

1.  Introduction 

In  some  form  Gauss  elimination,  or  GE  as  we  shall  often  abbreviate  it  here,  is  probably 
the  most  widely  used  single  computational  tool  in  scientific  computing  for  a  very  wide 
range  of  underl}dng  problems.  These  include  solution  of  partial  differential  equations, 
linear  programming  and  least  squares  approximation  in  its  various  guises  such  as  signal 
processing  and  linear  regression.  The  basic  problems  for  which  it  is  used  are  the  solution 
of  linear  systems  of  equations  (including  obtaining  the  solution  space  for  underdetermined 
svstems) ,  computation  of  matrix  determinants,  determination  of  matrix  rank  or  detection 
of  singularity. 

This  report  is  intended  to  provide  a  unified  treatment  of  one  of  the  most  important 
tools  of  linear  algebra  and  of  scientific  computing.  Gauss  elimination  is  usually  to  be  found 
as  a  major  topic  in  texts  on  linear  algebra  (where  the  emphasis  is  on  its  theoretical  basis), 
numerical  analysis  (with  an  emphasis  on  its  practical  implementation,  paying  attention 
to  questions  of  roundoff  error  and  numerical  stability),  parallel  and  vector  computing  (as 
an  important  illustration  of  the  potential  power  of  parallel  computers).  See  [l],  [3],  [7], 
[17],  [18]  for  examples. 

Another  aspect  which  is  usually  ignored  in  those  particular  texts  is  the  interplay 
between  GE  and  the  computer  arithmetic  being  used.  Roundoff  error  analysis  covers 
traditional  floating-point  computation  for  real  (or  complex)  systems.  For  integer  systems, 
however,  the  problem  is  not  roundoff  error  but  the  growth  in  the  matrix  elements  as  the 
elimination  proceeds.  Recently  this  aspect  must  be  extended  beyond  computer  arithmetic 
to  include  questions  related  to  the  computer  algebra  system.  In  addition  to  the  texts 
mentioned  above,  GE  is  included  as  a  symbolic  routine  in  computer  algebra  systems  such 
as  Maple^  and  Mathematical.  In  this  context,  GE  is  not  necessarily  restricted  to  matrices 
wdth  numerical  entries,  but  includes  matrices  with  entries  in  other  algebraic  rings  such  as 
algebraic  polynomials.  See  [2],  [27]  for  example.  (Also  note  Maple  linear  algebra  library 
functions  such  as  pffge  which  performs  “fraction-free’  GE  over  polynomial  rings.) 

^  Maple  is  a  registered  trademark  of  Maple  Waterloo  Software,  Inc. 

^Mathematica  is  a  registered  trademark  of  Wolfram  Researdi,  Inc. 
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In  this  report  we  begin  with  a  description  of  the  basic  GE  algorithm  and  some  vari¬ 
ations  of  this.  Section  2  also  contains  a  brief  review  of  the  main  applications  and  of  the 
complexity  and  error  analyses  for  these.  Section  3  is  concerned  with  integer  GE.  The  re¬ 
quirements  of  integer  arithmetic  raises  some  different  issues.  Roundoff  error  is  no  longer  a 
problem,  of  course.  Since  the  integers  are  not  closed  under  division,  we  describe  first  the 
division-free  form  of  the  al^rithm  and  then  the  questions  raised  by  this  -most  notably, 
the  growth  of  the  dynamic  range  that  is  needed  and  what  if  any  pivoting  strategy  should 
be  used  in  an  integer-arithmetic  setting.  Some  of  these  issues  were  previously  discu^ 
for  the  specific  context  of  residue  number  systems  in  the  series  of  reports  and  papers  [10]  ^ 


[ll]  [21]  [22],  [23]. 

Sections  4  and  5  deal  with  GE  for  matrices  with  entries  in  other  fields  or  rings.  Specif¬ 
ically,  Section  4  discusses  the  use  of  rational  arithmetic  while  Section  5  is  concerned  with 
more  general  rings  such  as  rings  of  polynomials.  The  question  of  solvability  is  addressed. 

Section  6  is  devoted  to  the  impact  of  various  forms  of  parallel  computer  archit^tures 
on  the  implementation  and  efficiency  of  GE.  Vector  and  array  architectures  are  considered. 
The  use  of  such  systems  for  performing  veiy^  long  integer  arithmetic  for  standard  integer 
GE  is  also  considered  here.  In  Section  7,  we  consider  GE  as  a  tool  for  a  specific  practical 
problem  -  namely  the  determination  of  the  effective  rank  of  a  matrix  whose  entries  are 
contaminated  by  noise.  One  of  the  key  aspects  here  is  setting  the  appropriate  tolerance 
level.  This  question  is  discussed  in  some  detail. 


1.1.  Problem  statement  and  notation.  Except  where  specificaUy  stated  we  shall 
consider  a  square  n  x  n  matrix  A  although  some  of  the  problems  and  solutions  are  sim¬ 
ilarly  valid  for  rectangular  matrices.  Elements  of  the  matrix  A  will  be  denoted  by  aij. 
Its  determinant  is  denoted  det>l  and  its  rank  by  r{^).  The  three  basic  problems  we 
are  concerned  with  here  are  the  solution  of  systems  of  equations,  computing  detj4  and 
determining  r  (>l).  For  definitions  of  any  of  these,  see  a  standard  text  such  as  [l]. 

In  the  case  of  systems  of  equations,  we  use  the  notation 


Ax  —  h 


and  the  elements  of  the  unknown  and  right-hand  side  vectors  are  denoted  Xi,  bi. 

The  matrices  L  =  {kj)  and  U  =  denote  lower  and  upper  triangular  matrices 
respectively.  Typically  they  will  be  lower  and  upper  triangular  factors  of  >1  so  that 
A  =  LU.  In  cases  where  any  pivoting  has  been  used  L,  U  are  the  factors  of  a  permut^ 
version  of  A.  That  \s  LU  -  PA  where  P  is  a  permutation  matrix  indicating  the  order  in 
which  the  rows  of  A  have  been  used.  (In  a  practical  implementation  of  GE,  the  matrix  P 
would  usually  be  stored  in  the  form  of  a  permutation  or  pivot  vector.  We  use  the  notation 
P  for  either  form.) 

In  the  context  of  either  error  analysis  or  matrices  contaminated  with  noise  we  shall 
use  E  for  the  error  or  noise  matrix.  In  case  we  wish  the  noise  level  to  be  explicit,  we  use 
Eff  to  denote  a  matrix  whose  elements  are  (independent,  identically  distributed)  normally 
distributed  random  variables  with  mean  0  and  standard  deviation  <7.  That  is  each  element 
of  Ea  is  drawn  from  N  (0,a). 

The  various  arithmetic  and  algebraic  systems  in  which  computation  is  taking  place 
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R  the  real  numbers 

Z  the  integers 

Q  the  rationals 

F  the  floating-point,  fl.p.,  number  system 

IL[x]  algebraic  polynomials  over  the  reals 

72.  a  general  algebraic  ring 

Clearly,  F  covers  a  wide  range  of  different  wordlen^hs,  precisions  and  implementations  of 
floating-point  arithmetic.  For  most  of  this  work  the  details  of  the  particular  fl.p.  system 
being  used  are  not  central  since  the  general  properties  are  similar.  Specific  fl.p.  systems 
will  be  detailed  as  needed. 


2.  The  basic  algorithms 

The  underlying  principle  of  GE  is  that  multiples  of  the  first  equation  (or  row  of  the 
matrix)  are  subtracted  from  all  subsequent  equations  (rows)  to  eliminate  the  first  unknown 
(element)  from  each  of  these.  The  process  is  then  repeated  to  eliminate  entries  below  the 
diagonal  of  each  column  in  turn.  The  system  of  equations  is  then  solved  by  substitution 
in  the  resulting  triangular  system.  Alternatively  the  determinant  of  the  original  matrix  is 
given  by  the  product  of  the  diagonal  entries  in  the  resulting  matrix.  The  rank  is  given  by 
counting  the  number  of  nonzero  elements  on  the  diagonal  of  the  final  array.  In  the  case 
of  solving  a  system  of  equations,  whatever  operations  are  performed  on  the  matrix  must 
also  be  performed  on  the  right-hand  side. 

Of  course,  these  statements  are  oversimplifications  of  the  true  situation  but  they  con¬ 
tain  the  essence  of  the  approach. 


Example  1.  We  illustrate  the  basic  idea  of  GE  /or  a  3  x  3  matrix. 


The  initial  matrix 


a  b  c 
d  t  f 
9  b  i 


is  modified  by  subtracting  d/a  times  the  first  row  from  the  second  and,  similarly,  g/a 
times  the  first  row  from  the  third  to  yield 


a 

b 

C 

a 

b 

c 

0 

e-ib 

= 

0 

e' 

f 

^  0 

h-H 

a 

i  — 

a  J 

0 

h' 

i’  _ 

say.  Next  h^/e'  times  the  second  row  of  the  modified  matrix  is  subtracted  from  the  third 


a 

b 

c 

a 

b 

c 

0 

e' 

/' 

= 

0 

e' 

r 

0 

p 

\ 

0 

0 

i"  _ 

from  which  we  obtain,  for  example, 

det  A 


Clearly  there  are  difficulties  in  the  event  that  either  of  a,  e'  are  zero.  If  we  had  been 
solving  a  system  of  equations  then  the  same  multipliers  would  have  been  used  to  adjust 
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the  right-hand  side  vector  so  that  the  original  system 


'  a  b  c 

’  X 

U 

d  e  f 

y 

— 

V 

9  h  i  ^ 

Z 

w 

would  be  reduced  to  the  triangular  system 


o 

1 — 

■  I  ■ 

“  1 

O 

O  O 

1 

_ 1 

1 - 

1 _ 

In  the  rest  of  this  section,  we  present  the  general  version  of  this  simple  GE  procedure 
and  discuss  some  of  the  difficulties  that  can  arise. 


2.1-  The  Ijk-form.  The  “ijk-form”  of  GE  is  just  the  general  n  x  n  version  of  the 
algorithm  used  in  Example  1  above. 


Algorithm  1  Basic  GE  ijk  form 

Input  n  X  n  matrix  A  (and  right-hand  side  b  if  solving  a  system) 

Compute 

for  t  =  1  to  n  —  1 

for  j  ^  +  1  to  r? 

771  OLjij  Q>ii 

Uji  :==  0 

hj  :=  hj  -  mhi  (if  solving  a  system) 
for  fe  =  r  -h  1  to  n 

Output  (modified)  matrix  A  (and  b  if  solving  a  system) 


Remark  1.  The  first  operation  in  tie  middle,  j-  loop  is  the  computation  of  the  appropri¬ 
ate  multiplier  so  that  the  apparent  repetition  of  the  division  d/a  m  Example  1  is  avoided. 
Clearly  the  algorithm  fails  if  an  =  0.  This  is  the  ext.reme  case  of  the  need  for  pivoting 
which  Is  discussed  below. 


Remark  2.  The  second  line  of  the  j-loop  sets  the  element  being  “zeroed”  directly  rather 
than  performing  the  arithmetic  which  should  yield  that  0. 


Pivoting.  Pivoting  is  the  name  given  to  altering  the  order  in  which  the  rows  are 
used  in  the  GE  algorithm.  The  simplest  cause  for  the  need  for  pivoting  is  the  occurren^ 
of  a  0  in  the  appropriate  diagonal  position  so  that  Algorithm  1  breaks  down.  In  t  e 
floating-point  environment  an  almost  equally  difficult  situation  is  created  by  a  vep’  smafl 
diagonal  entry  which  may  be  just  the  result  of  roundoff  error  in  a  quantity  which  would 
be  0  if  exact  arithmetic  were  being  used.  Such  a  pivot  element  can  result  in  large  errors 

in  the  computed  solution.  •  i  i 

The  usual  pivoting  strategy  is  partial  pivoting  in  which  the  curreirt  i-th  column  is 

searched  for  its  element  of  largest  magnitude  on  or  below  the  diagonal.  The  row  in  which 
this  occurs  is  then  interchanged  with  the  r-th  row  before  the  elimination  continues.  This 
form  of  GE  is  described  in  Algorithm  2. 
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Algorithm  2  GE  with  Partial  Pivoting  ijk  form 

Input  nxn  matrix  A  (and  right-hand  side  b  if  solving  a  system) 

Compute 

for  i  =  1  to  n  —  1 

find  p>i  such  that 

jopil  =  max  {|0ji|  :i<j<n} 

interchange  rows  i  and  p  (including  b  if  solving  a  system) 
for  j  =  t  +  1  to  71 


m  :=  Oji/aa 

aji  1=  0 
bj  :=  hj  -  mbi 
for  =  7  +  1  to  n 


(if  solving  a  system) 


Ojk  *—  Tnaj^. 

Output  (modified)  matrix  A  (and  b  if  solving  a  system) 

Remark  3.  In  practice  it  is  usually  simpler  to  keep  a  record  of  the  order  in  which  the  rows 
are  used  by  stoSng  a  permutation  vector  (which  is  equivalent  to  storing  the 
matrix  P)  If  the  interchange  is  performed,  it  is  of  course  only  necess^  to  interchange 
the  entries  for  columns  i  through  n  since  the  first  i  -  1  elements  of  both  rows  are  already 

zero. 

LU  factorization.  For  many  purposes,  it  is  desirable  to  make  better  use  of  the  work 
involved  in  GE.  By  storing  the  multipliers  used,  we  obtain  the  LU  factorization  of  the 
original  matrix  A.  That  is  we  find  a  lower  triangular  factor  L  and  an  upper  triangular 

one  U  such  that 

A  =  LU 

or,  in  the  case  of  pivoting,  L,  U  are  factors  of  a  permuted  version  of  the  matrix  A: 

PA  =  LU 

Algorithm  3  Basic  LU  Factorization  ijk  form 

Input  nxn  matrix  A 
Compute 

for  7  =  1  to  n  —  1 

for  j  =  7  -h  1  to  71 
aji  :=  Ojijoii 
for  /c  =  7  +  1  to  n 

CLjk  ■=  ^jk  ““  Q'ji^ik 
Output  (modified)  matrix  A 

Remark  4.  For  this  form  of  the  algorithm,  the  factors  L  and  U  are  stored  in  the  saine 
locations  as  the  original  matrix  A.  The  lower  factor  L  has  unit  diagonal  entries  and  so 
these  need  not  be  stored  explicitly. 

Remark  5.  This  choice  o/ a  unit  lower  triangular  factor  is  com'enieut  but  the  two  fac¬ 
tors  can  be  scaled  in  a  variety  of  ways.  This  choice  is  known  as  Doolittle  factoiization. 
Modifying  the  algorithm  to  yield  a  unit  upper  triangular  factor  is  the  Crout  reduction. 

Remark  6.  If  A  is  symmetric  the  factors  can  be  scaled  to  shai-e  the  same  diagonal  entries 
-  tJlis  is  then  called  the  CJjoIesky  factorization. 
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Reasons  and  advantages.  The  principal  advantages  of  the  LU  factorization  lie 
in  the  fact  that  if  multiple  systems  are  to  be  solved  with  the  same  coefficient  matrix 
then  the  factorization  can  be  done  just  once  and  each  system  solved  using  forivard  and 
back  substitution  loops.  This  is  then  much  cheaper  computationally  than,  for  example, 
inverting  the  matrbc.  This  is  detailed  in  Section  2.4  below. 

On  a  serial  computer,  inverting  a  matrix  (of  dimension  greater  than  2  x  2)  is  almost 

never  the  right  approach  to  a  problem. 

Pivoting.  The  LU  factorization  has  the  same  need  for  pivoting  as  does  the  basic  GE 
algorithm.  The  determination  of  the  appropriate  pivot  element  is  just  the  same  as  for  GE. 
It  is  necessary  to  keep  a  record  of  the  interchanges  in  order  to  make  the  corresponding 
changes  to  a  right-hand  side  vector  or  to  obtain  the  correct  sign  for  the  determinant. 

For  these  reasons  it  is  more  convenient  to  simply  store  the  permutation  matrix  (or 
vector)  P  for  subsequent  use  in  the  solution  process.  This  is  the  form  described  in  the 
following  algorithm. 

Algorithm  4  LU  Factorization  with  Partial  Pivoting  ijk  form 


Input  nxn  matrix  A 
Initialize  permutation  vector 
for  i  =  1  to  n 

P[»:]  :=i 

Compute 

for  7  =  1  to  n  —  1 

find  p>i  such  that 

{|ap[j)il  :i<j<n} 
interchange  P  [i]  and  P\p] 
for  j  =  ?■  -h  1  to  n 

apij]i  :=  apij]i/apii^i 

for  ^  f  +  1  to  n 
Output  (modified)  matrbc  A 


Remark  7.  The  *'IoweP'  /actor  has  elements  Uj  =  ^P\i]j  ^  '>  3-  The  "upper  factor 
has  Uij  =  apii]j  for  i  <j.  The  product  of  these  two  is  PA  where  P  is  now  regarded  as 
the  permutation  matrix  which  has  O^s  everywhere  except  for  l^s  in  the  positions 
or,  equivalently,  pij  =  Sp\i]j- 


Example  2.  Consider  the  4x4  matrix 


12  3  4 

2  2  3  4 

3  3  3  4 

4  4  4  4 


The  initial  permutation  vector  is  P  =  (L  2,3,4).  Clearly  the  pivot  element  in  the  first 
column  is  041  =  4.  So  P  is  now  (4. 2. 3.1)  and  the  matrix  resulting  from  the  first  step  of 

[1/4123^ 

1/2012 
3/4  0  0  1 
4  4  4  4 


the  elimination  process  is 


.  None  of  the  subsequent  elimination  steps 
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alter  this  matrix  -  but  they  do  change  the  permutation  vector  to  (4, 1,3,2)  and  finally  to 


Then,  we  have  L  = 


1 

1/4 

1/2 

3/4 


4 

2 

2 

3 


4 

4 

4 

4 


which  is  also  PA  = 


1 

'4444' 

II 

1  2  3 

1  2 

1 

%  .  »  ..  .  n 

0 

0 

0 

1 


1 

2 

3 

4 


and  then  LU  = 


4 

4 

4 


2.2.  Algorithmic  variations.  There  are  several  variations  on  the  basic  ijk-torm  of 
GE  which  have  merits  for  different  computing  environments.  We  concentrate  here  on 
just  two  of  these  variations.  For  simplicity  we  take  no  account  of  pivoting  h^e,  and 
we  restrict  to  Gauss  elimination  rather  than  the  factorization  of  the  matrix  j4.  (  e 
modifications  for  these  are  straightforward.)  The  rjfc-form  in  Algorithm  1  is  well-suited 
to  a  conventional  serial  computer  architecture  with  matrices  stored  by  rows.  For  other 
storage  schemes  or  processor  architectures,  alternatives  may  be  preferred.  A  much  more 
detailed  discussion  of  these  aspects  is  included  in  [18].  (The  reader  should  be  aware  of  a 
difference  in  notation  between  this  report  and  [18]:  there  the  pivot  row  is  denoted  by  ^ 
and  the  current  elimination  row  by  i  so  that  our  ijk- form  is  denoted  the  /czj-form  in  [18J.) 


The  ikj-form  for  column  vector  operations.  The  algorithm  for  the  ikj-form  is 
as  follows. 


Algorithm  5  Basic  GE  ikj  form 

Input  n  X  n  matrix  A  (and  right-hand  side  b  if  solving  a  system) 

Compute 

for  f  =  1  to  T?  —  1 
for  j  =  7  -f  1  to  n 

m.j  :=  aji/aa:  aji  :=  0 
bj  :=  bj  -  mjbi  (if  solving  a  system) 
for  A;  =  7  +  1  to  n 
for  j  =  7  +  1  to  n 

CLjk  0>jk  mjQ>ik 

Output  (modified)  matrbc  A  (and  b  if  solving  a  system) 

Remark  8.  Sow  the  innermost  hop  is  performing  a  “column  +  scalarx  column  opera¬ 
tion  where  the  scalar  is  the  element  aa-  of  the  pivot  row  and  the  column  it  multiplies  is 
the  column  of  multipliers  mj  computed  by  the  initial  j-loop.  The  corresponding  operation 
in  the  ijk^fortJi  is  a  ‘Vow  -h  scBla^rxrow  operation. 

Both  of  these  are  examples  of  (in  the  terminolog.v  of  [7])  saxpy  operations.  The  ikj- 
form  is  better  for  matrices  that  are  stored  by  columns.  The  choice  of  optimal  algorithm 
is  therefore  system  dependent. 
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Scalar  product  forms.  For  processors  which  are  designed  for  rapid  execution  of  a 
multiply-accumulate  (MAC)  operation  such  as  s  :=  s  +  a  x  b  it  is  d^irable  to  organize  an 
algorithm  to  utilize  this.  The  basic  linear  algebra  operation  which  is  particularly  intense 
in  its  use  of  MAC  operations  is  the  formation  of  a  scalar  product.  Further  variations  of 
GE  are  therefore  designed  to  use  scalar  product  operations  as  much  as  possible. 

To  see  how  this  fits  into  the  general  pattern,  it  is  easiest  to  look  at  the  innermost 
operation  of  the  LU  factorization  Algorithm  3.  If  the  order  of  the  loops  is  alter^  so  that 
the  i-loop  is  the  innermost  one  then  the  operation  being  performed  would  ind^  be  a 
scalar  product.  The  details  are  more  complicated  and  again  the  reader  is  referred  to  [18J 
(Appendbc  1)  for  a  full  explanation. 

Algorithm  6  Basic  LU  factorization  jki  (scalar  product)  form 

Input  n  X  n  matrix  A 
Compute 

for  j  =  2  to  n 
for  /c  =  2  to  j 
dj.k-l  '= 
for  i  =  1  to  fc  —  1 

CLjk  0>jk  “■  ^ji<^ik 
for  /c  =  j  +  1  to  n 
for  i  =  1  to  j  —  1 

ajk  •=  O'jk  ”  0,ji0>ik 
Output  (modified)  matrix  A 

Remark  9.  As  in  the  case  of  Algorithm  3.  the  output  matrix  contains  the  Doolittle 
factors  of  the  original  matrix. 

Remark  10.  Although  the  loop  structure  of  Algorithm  6  is  more  complicated,  the  arith¬ 
metic  operation  count  is  unchanged.  The  k-loop  is  separated  into  two  parts  corresponding 
to  the  lower  and  upper  triangles  of  the  matrix.  The  innermost  loop  in  each  of  the  parts 
is  using  MAC  operdtions  to  update  matrix  entries. 


2.3.  Applications. 

Solving  linear  systems.  The  most  frequent  application  of  GE  and  LU  factorization 
is  to  the  solution  of  a  linear  system 

Ax  =  b  (1) 


Then  the  GE  algorithms,  Algorithms  1  and  2,  result  in  an  equivalent  upper  triangular 


svstem  ^ 

t/x  =  b 


(2) 


whose  solution  is  the  same  as  that  of  (1).  This  system  is  then  solved  using  back  substitu- 
tion. 


Algorithm  7  Back  substitution  —  Row  form 

Input  n  X  n  upper  triangular  matrix  U,  n- vector  b 

Initialize  solution  x  :=  0 

Compute 
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for  i  =  n  downto  1 
for  j  =  z  +  1  to  n 

hi  1“  hi  XlijXj 
Xi  bi/Uxi 

Output  solution  x 

Remark  11.  Note  that  b  has  been  used  (in  place  ofb)  to  denote  the  right  hand  side  of 
(2)  in  this  algorithm. 

Remark  12.  This  version  of  back  substitution  has  been  arranged  in  such  a  way  that  the 
outer  loop  results  in  obtaining  one  furtier  component  of  the  solution  each  time.  The  inner 
loop  of  this  algorithm  subtracts  the  scalar  product  of  the  i^th  row  of  U  and  the  current 
solution  vector  from  the  i-th  component  of  the  right  -hand  side.  It  is  thus  well-suited  to  a 
computational  environment  which  is  designed  for  efficient  multiply-accumulate  operations 
and  the  matrix  is  stored  by  rows. 

For  architectures  which  are  well-suited  to  “vector  plus  scalar  x  vector’  operations  — 
these  are  the  saxpy's  of  [7]  —  a  column-oriented  version  of  this  algorithm  is  probably  to 
be  preferred. 


Algorithm  8  Back  substitution  —  Column  form 

Input  n  X  n  upper  triangular  matrix  U ,  n-vector  b 

Initialize  solution  x  :=  0 

Compute 

for  j  =  n  downto  1 

Xj  :=  bj/ujj 
for  2  =  1  to  j  —  1 

hi  hx  —  UijXj 

Output  solution  x 


Remark  13.  Here  the  inner  loop  subtracts  Xj  x  the  j-th  column  of  U  from  the  entire 
right-hand  side  vector  as  soon  as  it  is  available. 

In  the  event  that  pivoting  has  been  achieved  by  using  a  pivot  vector  (or  matrix)  such 
as  in  Algorithms  3  and  4  then  that  pivot  information  must  be  carried  through  to  the  back 
substitution  phase.  Essentially,  this  means  replacing  all  row  labels  in  these  algorithms 
by  F[i)  or  P[j].  For  the  column-form  of  the  algorithm  this  may  result  in  some  loss  of 
efficiency  since  data  will  no  longer  be  accessed  in  a  natural  order.  For  this  reason  it  may 
be  worth  performing  the  row  interchanges  for  architectures  where  this  is  the  preferred 
algorithm.  These  leist  comments  apply  equally  to  the  LU  factorization  solution  of  the 

system  (1).  ■  r\  u 

For  the  LU  factorization,  Algorithms  3  and  4  result  in  (a  permuted  version  of)  the 


system 


LUx  =  h 


(3) 


For  simplicity  we  ignore  the  effect  of  pivoting  here.  Solving  (3)  is  equivalent  to  soKing 
Ly  =  b  and  then  Ux  =  y.  The  first  of  these  requires  a  forward  elimination  process  which 
is  (almost)  entirely  equi\’alent  to  the  back  substitution  used  for  the  second  stage.  The 
only  differences  are  that  the  loops  run  in  the  opposite  directions  and,  because  L  is  a  unit 
lower  triangular  matrix  the  division  by  the  diagonal  element  is  not  needed. 
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Remark  14.  We  can  noiv  see  the  advantage  of  LU  factorization  for  solvmg  wuittp  e 
systems  with  the  same  coefficient  matrix.  To  solve  a  second  or  subsequent  system,  would 
require  only  forward  and  back  substitution  phases  to  be  repeated  since  the 
would  be  identical.  For  the  elemental  form  of  GE,  the  elimination  stage  would  ^so  need 
to  be  repeated.  As  we  see  in  Section  2.4,  this  is  the  most  expensive  part  of  the  whole 

process. 

Determinant  evaluation.  Whether  we  use  GE  or  LU,  the  evaluation  of  detA  is 
extremely  simple  once  the  elimination  phase  is  complete.  (We  are  taking  no  account  of 

any  Questions  of  error  analysis  or  stability  at  this  stage.)  ,  ^  r  j-  i 

In  the  absence  of  pivoting  the  determinant  is  simply  the  product  of  the  diagonal 
elements  of  U.  We  describe  the  algorithm  in  terms  of  the  LU  factorization. 

Algorithm  9  Determinant  evaluation  -  without  pivoting 

Input  nxn  matrix  A  r  .  •  i  • 

Compute  LU  factorization  by  Algorithm  3  (or  any  of  its  equivalent  varia¬ 
tions) 

initialize  D  :=  1 
for  i  =  1  to  71 
D  D  ^  Q>ii 
Output  det  A  is  D. 

The  only  additional  difficulty  if  pivoting  is  used  is  keeping  track  of  the  interchanges 
that  have  been  made  (whether  explicitly  or  implicitly)  in  order  to  get  the  correct  sign 
for  the  determinant.  It  suffices  to  maintain  a  variable  s,  say  in  the  pivoting  algorithm 
which  is  initialized  to  s  :=  1  and  is  multiplied  by  -1  whenever  7  r  P  (m  Algorithm  2) 
or  P\i]^  P  \p]  in  Algorithm  4.  In  the  latter  case  the  multiplication  loop  in  Algorithm  9 

would  also  be  replaced  by 

for  7  =  1  to  71 

Z?  :=  D  *  Op[t]i 

Rank  and  singularity  detection.  Again,  neglecting  any  problems  created  by 
roundoff  errors  or  ill-conditioning  in  the  matrK,  once  either  GE  or  LU  (with  pivoting)  has 
been  completed,  singularity  is  detected  simply  by  seeking  a  0  entry  in  a  pivot  position. 
In  fact  it  is  sufficient  to  examine  o„n  since,  if  any  pivot  element  is  zero,  then  necessarily 

note  we  are  assuming  both  pivoting  and  exact  arithmetic  or  algebra  in  making 
this  statement.  In  practice,  for  floating-point  arithmetic  at  least,  this  is  not  sufficient  and 
more  care  must  be  taken  to  test  for  near- singularity. 

Extending  our  ideal  world  analysis,  by  counting  the  number  of  zero  pivots  obtain 
the  rank-deficiency  of  the  matrix.  Equi^•alently  the  number  of  nonzero  pivots  yields  the 
rank  of  A.  This  is  a  much  more  optimistic  claim  than  even  the  singularity  statement. 
Within  a  computer  algebra  system  or  with  exact  arithmetic  such  statements  are  correct. 
In  our  ideal  world,  the  rank  algorithm  is  simple. 

Algorithm  10  Rank  detection  —  for  exact  arithmetic  or  algebra  systems 

Input  n  X  n  matrix  A  ... 

Compute  LU  factorization  with  pivoting  by  Algorithm  4  (or  any  of  its  equiv- 

alent  variations) 
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initialize  r  :=  0 
for  i  =  1  to  71 

if  api<)t  r  0  then  r  :=  r  +  1 
Output  rank  A  is  r. 

In  the  domain  of  real  numerical  computation,  the  question  of  rank  determination  is 
more  difficult.  The  effect  of  roundoff  error  on  GE  has  been 

of  that  analysis  is  summarized  in  Section  2.5.  Packages  such  as  MATLAB  (see  [16],  for 
example)  include  rank  as  a  function  in  their  computational  library  and  u^  a  c^efully 
computed  tolerance  dependent  on  parameters  of  the  computer  system  to  det^ine  the 
“true”  rank  of  the  original  matrbc  from  its  Singular  Value  Decomposition  (SVD)  [7].  In 
the  event  that  the  matrbc  entries  are  themselves  subject  to  error  —  perhaps  experimental 
error  or  the  effect  of  noisy  data,  for  instance  -  the  problem  becomes  more  difficult  to 
handle  reliably.  This  is  discussed  further  in  Section  2.5  and  then  in  Section  7. 

The  solution  space  for  underdetermined  systems.  In  the  case  of  underdeter¬ 
mined  systems  of  equations,  whether  this  is  the  result  of  having  fewer  equations  than 
unknowns  or  of  rank  deficiency  in  the  matrix,  a  set  of  vectors  spanning  the  solution  space 
is  easy  to  obtain  from  the  LXJ  factorization  of  the  matrix.  ^ 

The  backward  substitution  must  be  performed  the  same  number  of  times  as  the  rank 
deficiency.  Thus  if  the  ti.  x  n  matrix  has  rank  r,  we  solve  the  smaller  sj’stem  n  —  r  times 
each  time  with  r  of  the  unknowns  specified.  To  guarantee  linearly  independent  solutions, 
it  suffices  to  use  the  standard  basis  vectors  ei,  62, . .  • ,  e„_r  in  turn  to  specify  the  values  of 
Xr+i,a:r+2,..-,x„.  That  is,  we  first  set  . .  ,Xn)  =  ^  q": 

the  remaining  unknowns.  The  process  is  repeated  for  (Xr+i ,  Xr-r2, . .  •  r  ^ 

and  so  on  until  we  have  used  (Xr+i,z:r+2,- •  •  =  (Oi- •• !)•  _ 

The  algorithms  are  simple  extensions  of  those  outlined  earlier.  Similar  basis  vectors 
could  be  used  in  more  abstract  settings  with  symbolic  computer  algebra  systems. 


2.4.  Elementary  complexity  analysis.  The  computational  complexity  of  floating¬ 
point  algorithms  is  often  measured  by  the  number  of  floating-point  arithmetic  operations 
that  are  required.  Ti-aditionally,  the  relative  efficiency  of  algorithms  was  measured  by 
counting  multiplications  and  divisions  and  neglecting  addition  and  subtraction  opera¬ 
tions.  For  modern  processors  the  time  required  for  floating-point  multiplication  is  not 
much  greater  than  that  for  addition  and  so  it  is  sensible  to  consider  all  arithmetic  opera¬ 
tions.  Division  still  costs  substantially  more  and  any  elementary  function  evaluations  are 
typically  yet  more  expensive.  The  floating-point  operation  counts  for  GE  are  well-known. 
They  can  be  found  in  almost  any  standard  text  on  numerical  analysis  or  computational 
linear  algebra  such  as  [3]  or  [7]. 


Operation  counts. 


TABLE  1  Floating-point  operation  counts  for  GE  solution  of  an  n  x  n  linear  system 
Ax  =  b. 


Operation 

Forward  elimination 

Back  substitution 

TOTAL 

+  or  — 

X 

/ 

It?  11 

In  -  1) 

i7i(n-  1) 

|71.{71.-  1) 

n 

•in  (n  —  1)  (2n.  +  5) 
|n(n-l)  (2n  +  .5) 
|n  (n  -1- 1) 

3  MATLAB  is  a  registerwl  tratlemark  of  Tlie  MatbWorks,  Inc 
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If  multiple  systems  with  the  same  coefficient  matrix  are  to  be  solved  then  the  overall 
operation  count  for  GE  is  obtained  by  multiplying  these  totals  by  the  number  of  systems. 
In  particular  the  inversion  of  the  matrix  A  would  therefore  result  in  multiplying  all  these 
totals  by  n  so  that  matrbc  inversion  by  GE  is  an  0(n'‘)  operation.  .  . 

The  corresponding  table  of  operation  counts,  Table  2,  for  Doolittle  factorization  mak« 
it  easy  to  see  the  practical  advantage  of  using  the  LU  factorization  whenever  multiple 
s^’stems  are  to  be  solved. 


TABLE  2  Floating-point  operation  counts  for  solution  of  j4x  —  b  using  Doolittle  LU 
factorization. 


Operation 

Factorization 

Forward  sub 

Back  sub 

TOTAL 

+  or  - 

X 

/ 

|n(n—  1)  (2n  —  1) 
|n(n—  l)(2n-  1) 
\n{n-  1) 

4n(n-l) 

|n(77  -  1) 

0 

ln{n  -  1) 
\nln-  1) 

n 

in  (n  -  1)  (2n  4-  5) 
|n(n-l)(2n-l-5) 
In  (n  4- 1) 

Now  in  order  to  solve  several  systems,  the  factors  of  the  original  matrix  are  already 
known  and  so  only  the  operation  counts  for  the  forward  and  backward  substitution  phases 
need  to  be  repeated. 


TABLE  3  Floating-p>oint  operation  counts  for  solution  of  m  systems  using  LU  factor¬ 
ization  and  GE 


Op 

Factorization 

1  system 

TOTAL  for  LU 

TOTAL  for  GE 

+  or  - 

X 

/ 

in  (n  —  1)  (2n  —  1) 
ln(n-l)(2n-l) 
in(n-  1) 

n  (n  ~  1) 
n(n  —  1) 

n 

^n(n-l)(2n-l) 
H-7nn(n—  1) 
iTi(n-l)(2Ti-l) 
-LTnn(n- 1) 

^n(7?.  —  1)  +  mn 

f  n  (n  -  1)  (2n  4-  5) 
^n(n  —  1)  (2n  4-5) 
^nin  +  1) 

It  is  immediately  apparent  that  the  savings  here  are  piotentially  great  for  large  systems 
with  multiple  right-hand  sides.  For  inversion  of  a  large  matrbc,  GE  requires  approximately 
2n^/3  flops  compared  to  8n^/3  for  LU  factorization. 


2.5.  Elementary  error  analysis.  In  this  section,  we  summarize  briefly  the  well- 
known  error-analysis  of  GE  and  LU  factorization.  For  a  more  complete  description  see  [3] 
or  [7].  For  all  the  usual  error  bounds,  we  need  the  notation  associated  with  the  standard 
vector  and  matrix  norms  and  the  condition  number  of  a  matrix. 

Condition  numbers  and  norms.  There  are  three  impoi-tant  vector  norms  that 
are  commonly  used  in  anlaj-sing  numerical  methods  for  matrbc  problems.  These  are  the 
L1.L2  and  Loo  norms  which  are  defined  for  x  €  R"  by 


llxlli  = 

H  l^«l 

(-») 

i=l 

l|x||2=  . 

(5) 

\ 

t=l 

ilxiloo  =  kil 

(6) 
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which  are  often  known  as  the  taxicab.  Euclidean  and  maximum  norms  respective!}. 
Any  vector  norm  has  an  associated  matrix  norm  given  by 


|1A||  =  max{||Ax||/l|x||} 

In  particular  the  Li  and  Loo  norms  can  be  obtained  by 

||A||i  =  m^^lOy| 

X 

(7) 

||Al|oo  =  mp^iavl 

(8) 

} 


which  are  respectively  the  maximum  column  (absolute)  sum  and  maximum  row  (absolute) 

sum.  .  ■  /  /i\ 

The  condition  number  of  a  matrbc  (relative  to  whichever  norm  is  being  used)  is  k(A)  — 
||j4j|  •  ||j4~^||.  It  is  used  as  a  measure  of  the  inherent  difficulty  of  a  linear  algebra  problem. 
The  larger  the  condition  number  the  more  difficult  it  is  to  get  a  numerically  stable  solution. 
Typically,  a  large  condition  number  implies  that  A  has  at  least  one  verj’  large  or  at  least 
one  very  small  eigenvalue.  This  will  often  result  in  large  accumulated  roundoff  errors 
in  the  solution  of  a  system.  However,  if  the  solution  to  such  a  system  has  a  very  small 
component  in  the  direction  of  the  eigenvector  associated  with  such  an  eigenvalue,  the 
solution  may  still  be  obtained  to  high  accuracy  (so  that  the  problem  is  weU-conditioned) 
despite  the  ill-conditioned  nature  of  the  matrix. 

Computing  the  condition  number  of  a  matrbc  is  of  comparable  difficulty  (and  compu¬ 
tational  complexity)  to  inverting  the  matrbc  in  the  first  place.  However  there  are  good 
algorithms  for  estimating  the  condition  number  which  may  be  used  to  estimate  the  number 
of  correct  digits  in  the  solution  to  a  linear  system.  These  are  discussed  in  [7]. 

Error  estimation.  The  conventional  way  of  measuring  the  error  in  the  solution  to  a 
linear  system  is  by  what  we  might  term  a  “relative  norm  error  ".  Recently  there  has  been 
extensive  work  on  reanalysing  methods  in  terms  of  a  componentwise  relative  error  which 
overcomes  some  of  the  overestimation  of  errors  which  may  arise  out  of  a  well-(X)nditioned 
problem  which  has  an  ill-conditioned  matrix.  (For  more  detail  of  this  approa,ch  seej4].) 

Suppose  that  the  exact  solution  of  Ax  =  b  is  x’  while  the  computed  solution  is  x.  We 
denote  by  r  the  residual  vector 

r  =  b  -  Ax 


The  relative  error  in  the  solution  can  now  be  bounded  by 


1  IWI  ^l|x-x-|| 
^{A)  l|b|l  -  l|x*ll 


<k(  41-!^ 
-  ^  Mlbll 


(9) 


(This  residual  vector  is  also  useful  for  carrying  out  iterative  refinement  of  the  computed 
solution  by  solving  the  system  Ay  =  r  using  the  LU  factors  and  then  adding  y  to  the 

original  solution  x.)  .  r  a-  r  ■ 

The  other  error  bounds  which  are  useful  here  give  estimates  of  the  effect  of  errors  in 

the  original  matrix  and/  or  right-hand  side  vector.  The  effect  of  an  error  in  b  is  estimated 
in  a  similar  manner  to  the  use  of  the  residual  above.  If  5b  is  the  error  vector  (meaning 
that  we  have  solved  Ax  =  (b  -I-  5b)  instead  of  Ax’  =  b)  then  we  obtain 


_LJM< 

k{A)  ||b||  - 


l|x-x’ 

llx’ll 


<  <A) 


IM 

l|b|| 
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which  we  see  has  just  the  same  form  as  (9).  errors  is 

The  corresponding  bound  for  the  situation  where  the  matrix  A  has  errors  is 


llx-x'l 


<  <A)'- 


These  various  bounds  can  be  combined  to  estimate  the  accuracy  of  a  final  solution  in  the 
“r„roTlndofr  erro^  .nd,  p»h,ps,  „oi«.  Again  aither  13]  [T]  to,  mce  detaila. 

SSSabk  im^STnce  put  on  the  detection  of  the  :^J“®Ynrifappli^t1rn?  T 

a  statistiL  analysis  of  the  SVD  (singular  value  decomposition)  alpnthm  to  obtain 
appropriate  tolerLce  levels  in  order  to  get  reliably  accurate  rank  information.  Mwe 
rSnth'  Gleeson  [6]  has  undertaken  a  study  of  the  viability  of  usmg  a  (precomputed) 
statistical  studv  of  coefficients  of  the  characteristic  polynomial  of  the  matrix  in  order  to 
identify  the  “signature”  of  a  matrix  of  given  rank.  The  idea  behind  these  studi^  k  to  try 
to  obt^n  the  rfnk  without  the  need  for  a  full  eigenvalue  or  singular  value  analysis  of  the 

Gauss  elimination  offers  a  cheap  alternative  to  these  apprc«ches  whk^  "f  M 

of  the  time.  (This  statement  is  contrary  to  that  made  in  [12]  which  was  ^  on  a  ^.^^eful 
and  LTericallv  unstable  implementation  of  GE.  However  in  the  case  of  an  ^^l-^ondition^ 
m^The  GE  approach  still  fails.  Briefly  the  use  of  GE  for  effective  rank  determination 
is  based  on  the  oServation  that  the  noisy  matrix  is  almost  surely  nonsin^lar  as  a  r^ult 
of  the  combination  of  roundoff  and  noise.  Therefore  all  entries  on  the  diagonal  of  the 
upptr  triangular  factor  will  be  nonzero  and  the  algorithm  will  be  successfully 
Recounting  the  number  of  these  diagonal  entries  which  are  greater  than  some  well-chosen 

threshold  value,  the  “correct”  rank  is  obtained.  ^  o  i.  ^ 

The  “correct”  tolerance  for  this  algorithm  appears  to  be  approximately  nS  where  n 
is  the  dimension  of  the  matrix  and  S  is  a  bound  for  the  elements  of  the  noise  matrix^ 
With  this  tolerance  most  matrix  rank  problems  are  correctly  resolved  -  but  most 

nec^rfi^^fficiei  has  a  built-in  rank  function  which  can  take  a  tolerance  param¬ 

eter  ThT^goriVhm  is  based  on  the  singular  value  decomposition  of  the  matrix  and  is 
to^refore  mucri^^  expensive  to  compute.  There  remains  a  question  of  what  toleraiice 
to  ui  altruEh  again  nS  appears  to  be  a  good  choice.  This  algorithm  is  rehabk 

provided  that  the  tolerance  is  not  comparable  with  the  entries  (and  singular  of  A 

An  alternative  algorithm  [24]  based  on  the  use  of  least  squares  approximations  of  rows  of 
4  bv  linear  combinations  of  other  rows  is  currently  under  iny«tigation  in  anothei  NAW  C 
project  and  appears  to  have  some  promise  as  being  both  reliable  and  of  comparable  cost 

to  GE. 

3  Gauss  eli.mination  over  the  integers 
For  computation  over  the  integers,  it  is  necessary  to  make  changes  to  the  basic  GE  alg^ 
rithms  dL:ribed  in  Section  2.  The  most  obvious  cause  is  the  fact  that  the  integeis  aie  n^ 
closed  under  division.  This  has  the  side  effect  that  the  magnitudes  of  inte^rs  general^ 
durhfg  the  computation  can  gl  ow  rapidly.  The  range  of  integer  values  required  or  a^allable 
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is  known  as  the  dyru^mic  range.  Overflowing  the  integer  range  in  bin^> 
often  results  in  the  phenomenon  known  as  “integer  wraparound  (s^  P).  Chapter  1  lor 
example)  which  is  the  effect  of  the  “clock"  arithmetic  modulo  2  where  A  is  the  integer 
wordlength  in  bits.  Further  complications  that  arise  out  of  integer  arithmetic  (however 
is  oerformed)  include  choosing  a  good  pivoting  strateg>-. 

^In  a  Computer  Algebra  System  (CAS)  such  as  Maple  [2]  or  j 

these  problems  are  avoided  at  the  expense  of  computational  speed  since  exart  a>^^hmetK 
with  verj^  long  integers  can  be  performed  in  software.  Howev«  such 
very  slow  when  the  arithmetic  wordlength  gets  long  and  the  rate  o  gr 

and  arithmetic  formats  for  the  integers.  The  r^idue  number  si^ms  RNS  are  well-^i^ 
to  some  of  these  tasks.  In  such  a  sj'stem,  an  integer  is  represented  by  its  ^dulo  a 

number  of  different  prime  numbers.  The  advantage  here  is  that  dynamic 

range  is  achieved  by  extending  the  set  of  basis  primes  being  used.  ^ 

a  nLral  short  wordlength  parallelism  which  avoids  the  slowdown  caus^  by  very  long 
wordlength  arithmetic.  However  it  brings  with  it  other  difficulties  ^hich  are  less  easily 
resolved  Even  if  one  integer  divides  another  and  both  are  within  the  ^namic  range  of 
the  system,  there  is  no  simple  division  algorithm  which  returns  this  result;  range  checking 
is  (at  best)  very  difficult;  comparison  is  not  an  RNS  operation  vi,  +• 

In  this  section,  we  discuss  the  implementation  and  use  of  GE  using  inte^  arithmetic 
Our  primary  focus  will  be  on  binary  integer  arithmetic.  The  use  of  RNS  anthme^c 
within  the  specific  context  of  adaptive  beamforming  was  discussed  in  a  seri«  of  reports 
and  papers  [9],  [10],  [11],  [21],  [22],  [23].  Specific  reference  to  using  RNS  arithmetic  will 
be  included  only  where  it  is  helpful  to  understanding  in  the  present  context. 

3.1.  Division-free  Gauss  elimination.  The  simplest  way  of  modifying  GE  to  inte¬ 
ger  arithmetic  is  to  eliminate  the  divisions  by  performing  “cross-multiplications  '  between 
the  rows  of  the  matrix.  This  is  the  form  which  generates  the  greatest  rate  of  ^oirth  in 

the  dynamic  range.  This  corresponds  to  the  transformation  of  2  x  2  matrices  ^  ^ 

This  is  achieved  by  multiplying  the  second 


instead  of 


a  h 
0  ad -- he  ^ 
row  by  a,  and  subtracting 


a  h 

,  0  d-  6(c/a)  I  ,  1 

rom  it  c  times  the  first  one.  The  overall  divisionless  GE  algo- 


row  DV  a,  ana  ^  — - •  r 

rithm  consists  of  repeating  this  for  all  the  appropriate  submatrices  as  in  the  basic  forms 

of  the  algorithm  in  Section  2. 

As  with  the  real  arithmetic  forms  of  GE  there  are  many  ways  of  arranging  the  o^ 
erations.  We  concentrate  here  on  just  the  simplest  -  the  fjfc-form  corresponding  to 

Algorithm  1. 

Algorithm  11  Division-free  GE  ijk  form 

In-put  n  X  n  integer  matrix  A  (and  right-hand  side  b  if  solving  a  system) 


CompxLte 


for  i  =  1  to  n  -  1 
for  j  =  r  +  1  to  n 

bj  :=  aabj  -  CjA  (if  solving  a  system) 
for  A:  =  t  -f  1  to  n 

Ojk  "■  ^ji^ik 
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Uji  :=  0 

Outvut  (modified)  matrix  A  (and  b  if  solving  a  system)  ,  .  ,  , , 

or  determinant  crienlation.  Three  dtorssod  later  m  Sf  ^ 

divieionlees  eUmination  on  the  resultine  matrix  elements  has  been  anal>5ed. 

3  2  Growth  to.  dynamic  range.  The  qinetion  of  the  range  grovrth  in 

Ge'wos  widreesed  in  some  detail  in  »  [11),  N  in  order  to  “ 

RNS  arithmetic  within  the  context  of  wiaptive  b^formmg.  ^  * 

snmmarizing  these  resoHs  general  int^ianthmet^In^erto^^^^^  totn^ry 
Sli^'thjMSx'^rte  misini  from  the  divisionless  algorithm  am^  corresponding 

remark  that  the  initial  range  does  not  necessarily  match  the  available  dynamic  range.  e 

symmetric  range  simplifies  the  analysis  of  the  powth. 

0)  0  .  ^  in  slu 

The  transformation  of  the  2x2  matrix  ^  ^  0  ad  — be  \ 

1  +  nW  A/*  c  r  M^l  If  M  =  2^  -  1,  this  implies  that  the  dynamic  range 

has'inLtJ^ bom  a  U'r^  minimum'wordlength  of  ft  +  1  biu  to  2K  +  2 
bto.  mis  same  exponential  rate  of  gttnvth  is  possible  at  every  stage  of  the  enter  loop 

^'®Th»”h“initial  required  wordlength  is  (approximately)  doubled  N  -  1  times  during 
thc  ^Xi^i^  GE^mination  lor  a%  x  hf  matrix.  The  final  wordtogth  needrf  » 
therefore  around  2"-'  (ff  + 1)  I‘  ^  oast-  to  see  that  this  would  very  quickly  exhaust  any 

normally  available  dynamic  range.  .  f _ 7  7I 

For  example,  if  were  just  3  so  that  the  range  is  r«tri^  to  just  [_  .  7] 

V  =  6  the  final  dvnamic  range  would  need  a  wordlength  of  at  least  2  x  4  -  laJS 
b£  which  4  alreadv  double  the  length  of  any  commonly  found  b-lt-m  integer  format^ 
Ctorly  more  Realistic  range  for  the  initial  data  and  larger  matrix  dimensioto  this 

eroinh  would  very  quicklv  exceed  all  plausible  ranges  even  for  software  implementation. 

intetm  arithmetic  as  in  |10)  the 
rapid  than  the  above  analysis  suggests  and,  perhaps  more 

ZTth  is  achieved  in  realistic  examples.  Clearly  it  becomes  necessary  to  find  ways  of 
Sstrkting  this  growth.  We  discuss  this  aspect  shortly  but  first  need  to  get  a  clear  picture 
of  the  comparative  magnitudes  of  matrix  elements  resulting  from  the  divisionless  algorithm 

and  that  using  division. 

Comparison  between  matrix  entries  with  and  without  divisions.  In  this 
subsection  we  look  in  some  detail  at  the  relations  between  the  matrix  elements  gener¬ 
ated  bv  the  divisionless  form  of  GE  and  those  that  would  be  obtained  with  dn  isioiis. 
Thfs  i^iniportant  for  the  completion  of  the  various  applications  especially  ^ 

2;:rminant  which  was  specially  straightforward  for  the 

usine  division  In  order  to  make  this  comparison,  we  assume  that  the  dnisions  can  be 
performed  exactly  in  some  rational  arithmetic  system.  Since  these  divisions  onlj  occur 


NAWCADPAX-96-194-TR 


Gauss  Elimination: 


Workhorse  of  Linear  Algebra 


within  an  entirely  theoretical  version  of  the  algorithm,  the  comparisons  remain  valid  and 

In  [10]  it  is  stated  that  the  the  growth  of  the  matrix  elements  is  such  that  the  final 
value  of  aNfi  at  the  conclusion  of  the  elimination  phase  is  detA.  This  turns  out  to 
be  overoptimistic  in  general.  The  potential  grow’th  is  much  greater  than  this  although 
it  becomes  apparent  from  the  analysis  that  the  algorithm  can  be  modified  to  have  this 
result.  First  we  need  some  notation  to  distinguish  between  elements  resulting  from  the 
different  forms  of  GE  under  consideration.  For  simplicity  here  we  shall  completely  ignore 
any  pivoting  and  shall  assume  that  the  algorithm  does  not  break  down  due  to  a  zero  pivot. 

Following  convention,  we  denote  by  o,j  the  elements  of  the  original  matrix.  All  refer¬ 
ence  will  be  to  the  simplest  tjfe-form  of  GE  described  by  Algorithm  1.  By  the  i-th  stape  of 
the  algorithm  we  mean  the  outermost  loop  for  the  value  i  which  performs  elimination  be¬ 
low  the  diagonal  in  the  i-th  column.  The  entries  resulting  from  the  i-th  GE  with  divisions 
are  denoted  by  The  corresponding  entries  for  the  divisionless  form  Algorithm  11 
will  be  denoted  by  l»j‘2  (j,  k>i).  The  “final”  values  in  the  two  algorithms  are  therefore 

given  by  a%'^\  (j  =  0,1, -1-,  k>j)  where  =  ai*. 

To  see  the  relations  between  the  entries  snd  j 

d]^[o  od-6c]  [o  d-6(c/a)]- 

in  the  current  notation, 


=  ad  —  6c  =  a  [rf  —  6(c/a)]  =  aua^ 


and  that  the  corresponding  relation  would  hold  for  all  other  matrix  entries  resulting  from 
the  first  stage.  That  is 

=  i2<j,k<N)  (10) 

In  particular,  we  see  that  is  the  determinant  of  the  top-left  2x2  submatrix.  It  will 
be  convenient  to  denote  this  by  dj-  In  general,  we  shall  denote  by  d;  the  determinant  of 
the  top-left  i  xi  submatrix  of  A. 

The  next  stage  of  the  elimination,  is  essentially  the  same  as  this  first  stage  operating 
on  the  bottom-right  (A^ - 1)  x  (A" - 1)  square  submatrix  —  the  “active  matrix^’.  It  follows 
that  there  is  a  further  scaling  of  all  affected  entries  by  the  pivot  element  Thus,  we 
deduce  that  .  . 

(3  <j,k<  N)  (11) 

and  continuing  in  this  way  we  obtain  the  general  relation: 


~  n.11622  "‘b. 


{i<j,k<N) 


Each  “final”  divisionless  entry  is  its  corresponding  value  from  Algorithm  1  scaled  by 
the  product  of  all  the  final  diagonal  entries  of  the  divisionless  algorithm  above  it. 

Remark  15.  note  immediately  that  this  analysis  suggests  that  entries  in  the  active 
matrix  have  common  factors  which  have  been  included  by  the  algorithm  itself.  These 
represent  one  obvious  waj'  of  reducing  the  range  growth  and  altering  the  algorithm.  This 
modification,  which  is  presented  beloiv,  w'ill  result  in  each  diagonal  element  being  the 
determinant  di  of  the  principal  minor  of  the  appropriate  dimension. 
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3.3.  Reducing  the  range  growth.  Before  considering  how  to  take  full  advantage  ^f 
the  enforced  common  factors  which  result  from  the  divisionless  algorithm,  we  lo^  at  the 
simplest  range  reduction  technique  resulting  from  the  removal  of  unnecessary  factors. 

Greatest  common  divisors.  The  simplest  approach  is  just  to  divide  out  any  com¬ 
mon  factors  in  the  “cross-multiplication”  operations  which  are  at  the  h^  of  the  dm- 
sionless  GE  algorithm.  Again,  this  idea  is  easily  understood  by  considering  the  2  x 

^  Suppose  that  a,c  have  a  common  factor  p 


elimination  beginning  with  A  — 


so  that  there  exist  integers  a,  7  such  that  o  =  pa,  c  -  P7 
'a  b  '' 

0  ad  — be 


transformation  to 


Then  the  usual  divisionless 
can  be  replaced  by  subtracting  7  times  the  first  row 


from  a  times  the  second  to  yield 


since  ac  -  07  =  ap^  -  pa7  =  0. 


a  b 

0  ad-bj  j 

This  suggests  that  it  is  desirable  to  find  the  greatest  common  divisor,  gcd(a,c),  of  a,  c 
before  proceeding  with  the  elimination.  This  same  principle  can  be  applied  at  each  step 
of  the  process  and  it  can  be  applied  in  more  general  settings  than  just  integer  arithmetic. 
For  example,  this  idea  is  incorporated  into  the  polynomial  fraction-free  Gauss  elimination 

functions  for  Maple.  .  .  ,  ,  11  1  \  + 

Further  simplification  may  be  achievable  using  the  ged  in  the  (perhaps  unlikely)  event 

that  a  complete  row  of  the  matrbe  at  some  stage  has  a  nontrivial  common  factor.  This  too 
can  be  divided  out,  reducing  the  range  growth  in  subsequent  stages  of  the  elimination. 
Of  course  if  the  goal  is  to  obtain  det  A  then  a  record  of  (the  product  of)  these  factors 
must  be  kept  to  multiply  the  final  result.  Finding  the  ged  of  a  complete  row  of  a  nriatrix 
of  anything  more  than  very  small  dimension  may  be  at  a  price  which  does  not  justify  the 

The  simplest  technique  for  finding  gcd(a,  c)  is  the  Euclidean  algorithm  which  is  easily 
implemented  as  follows: 


Algorithm  12  Euclidean  algorithm  for  integer  ged 

Input  Positive  integers  c  <  a 
Initi-olize  q  ;=  c,  tmpc  :=  a 
Repeat 

tmpa  :=  tmpc;  tmpC  :=  q 
q  :=  tmpa  mod  tmpC 
until  q  :=  0 

Output  gcd(a,  c)  is  tmpc 

This  can  be  used  within  the  following  modified  divisionless  GE  algorithm. 

Fraction- free  algorithm.  A  fraction-free  version  of  GE  can  be  defined  using  the 
ged  algorithm  above  to  reduce  the  range  growth  effect.  The  resulting  algorithm  still  needs 

no  divisions. 

Algorithm  13  Fraction-free  GE  using  greatest  common  divisors  ijk  form 

Input  n  X  n  integer  matrix  A  (and  right-hand  side  b  if  solving  a  system) 

Comp  tile 

for  7=1  to  71—1 
for  j  =  7  -h  1  to  7? 
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p  :=  gcd(aii.aji) 
a  :=  Gu/p;  y  :=  aji/p 
bj  :=  abj  -  76,  (if  solving  a  system) 
for  &  =  t  +  1  to  n 

ajk  :=  Ckajk  -  7^*^ 
aji  :=  0 

Output  (modified)  matrix  A  (and  b  if  solving  a  system) 

Remark  16.  Since  we  iiave  no  advance  knowledge  of  the  existence  of  nontriviai  factors 
(that  isp>  1)  this  algorithm  does  nothing  to  moderate  the  worst  case  range  growth  of 
the  analysis  in  Section  3.2.  However,  if  the  wordlengths  grow  according  to  the  needs  of 
the  particular  application,  it  is  likely  that  some  savings  wUl  materialize. 

Using  the  gcd  of  complete  rows  can  be  incorporated  into  the  algorithm  at  anj’  stage. 
Potentially  such  factors  could  exist  in  any  row  of  the  active  matrix  and  such  factors  could 
be  sought  in  every  row  at  ever)’  stage  of  the  elimination.  The  algorithmic  changes  are 
straightforward  but  the  cost  is  likely  to  be  too  high  and  so  we  do  not  elucidate  further 
here  —  except  for  one  very  important  special  situation  which  does  arise  as  a  result  of  the 

elimination  process  itself. 

Prom  (11)  it  follows,  in  particular,  that 

l>3f  =  <4?  =  ®ii43 

where,  we  recall,  4  =  an42  43  ‘s  the  determinant  of  the  principal  3x3  minor.  We 
remark  that  da  is  an  integer.  It  follows  that  b^^  has  the  factor  On.  Using  this  same 
reasoning,  we  see  that  bff}  has  a  factor  On  for  every  j,k  ^  3  since  each  such  element  is 
On  times  the  (integer)  determinant  of  a  3  x  3  minor.  This  known  factor  can  easily  be 
removed  prior  to  the  next  stage  of  the  elimination. 

Similar  factors  are  introduced  into  the  active  matrix  at  each  subsequent  stage.  These 
too  can  be  divided  out.  The  factors  introduced  at  the  subsequent  stages  of  the  elimination 
are  the  (modified)  diagonal  entries.  Since  these  are  known  factors,  there  is  no  need  for 
any  calls  to  a  gcd  function. 

The  removal  of  these  factors  has  an  obvious  effect  on  the  range  growth  in  subsequent 
stages  of  the  elimination  even  though  the  worst  case  analysis  is  unchanged.  The  grow-th 
in  the  required  wordlength  would  need  to  be  computed  “on-the-fly '  in  order  to  take 
advantage  of  this  saving. 

The  division  by  these  factors  is  incorporated  into  Algorithm  14  below.  \\e  observe 
that  this  algorithm,  like  Algorithm  13,  is  fraction-free  but  not  division-free.  All  divisions 
that  are  performed  are  integer  operations  with  exact  integer  results. 

Algorithm  14  Fraction-free  GE  ijk  form 

Input  n  X  n  integer  matrix  A  (and  right-hand  side  b  if  solving  a  system) 
Compute 

for  T  =  1  to  n  —  1 
for  j  =  i  -h  1  to  n 
i}j 

for  /:  =  ?  -h  1  to  n 


(if  solving  a  system) 
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Ojfe  :=  Ciii<ijk  ~~  (>'ji^ik 

Oji  :=  0 

if  r  >  2  then  (removal  of  common  factors) 

for  j  =  i  +  1  to  n 

bj  :=  bj /ai-i,i-i  (if  solving  a  system) 

for  /c  =  i  + 1  to  Ti 

ajk  fl'jfc/O't-i.t-i 

Output  (modified)  matrix  A  (and  b  if  solving  a  sj’stem) 

To  gain  some  insight  into  the  saving  that  results  from  this  algorithm,  considCT  just  the 
final  value  of  a^N.  In  the  division-free  Algorithm  11,  the  analysis  of  Section  3.2  yields, 

using  (12)  for  i  =  N  —  I- 


h(l) 


which  in  turn  yields 


r.W 


(N-2)^  JN-l) 


"NN 


=  dN 


Oi1®22  '  "  ^N-1,N-1^NN 

.(A)!'"-" 


{4-=  [‘S 


/  iV-2  h(^-3)  \ 

<  Oji  1^622  j  ■  "  ®A--2,A’-2  J 


J 


jXN-Z)  \ 


(13) 


The  corresponding  final  value  for  Algorithm  14  is  just 

O’NN  ~  det  A  ~  dft’  (I^) 

which  is  typically  very  much  smaller  since  these  additional  factors  have  been  remov^ 
during  the  modified  elimination.  Indeed  it  turns  out  for  this  Algorithm  14  that  at  the 
conclusion  of  the  elimination  we  have 

aii  =  di.  (f^) 

Example  3.  Comparison  of  Algorithms  11  and  14  for  a  4x4  matrix.  The  results  of  exact 
arithmetic  in  Algorithm  1  are  also  included  for  comparison  and  illustration. 


Let 


A  = 


8  7  4  1 
4  6  7  3 
6  3  4  6 
4  5  8  2 


Algorithm  1  yields  the  following  sequence  of  modified  matrices: 


8  7  4  1 

0  5/2  5  5/2 

0  -9/4  1  21/4 

0  3/2  6  3/2 


8  7  4  1 

0  5/2  5  5/2 

0  0  11/2  15/2 

0  0  3  0 


[8 

7 

4 

1 

0  5/2 

5 

5/2 

0 

0 

11/2 

15/2 

0 

0 

0 

-45/11  _ 

45/11) = 

-450. 

The  division-free  Algorithm  11  gives: 


NAWCADPAX-96- 194-TR 


Gauss  Elimination:  Workhorse  of  Linear  Algebra 


[87411  [874  l]  874  1 

0  20  40  20  0  20  40  20  0  20  40  20 

0  -18  8  42  0  0  880  1200  0  0  880  1200 

0  12  48  12  0  0  480  0  .  ®  ®  ®  —576000 

which  shows  the  very  rapid  growth  which  is  possible  with  this  algorithm. 

The  fraction-free  Algorithm  14  also  gives 

'87411  [874  l' 

0  20  40  20  0  20  40  20 

0  -18  8  42  0  0  880  1200 

0  12  48  12  _  [  0  0  480  0 

but  the  active  part  of  this  matrfac  is  then  divided  by  the  common  factor  8  to  yield 

'  8  7  4  1  ' 

0  20  40  20 

0  0  no  150 

0  0  60  0  _ 

which  in  turn  gives 

■  8  7  4  1 

0  20  40  20 

0  0  no  150 

0  0  0  -9000 

The  active  matrbc  is  now  just  the  bottom-right  element  which  needs  to  be  divided  by  the 
“common  factor”  20  to  yield: 

■  8  7  4  1  ■ 

0  20  40  20 

0  0  no  150 

0  0  0  -450 

The  final  values  of  the  diagonal  entries  can  easily  be  seen  (by  comparison  with  the  partial 
products  of  the  diagonal  of  the  final  upper  triangle  generated  by  Algorithm  Ij  to  be  the 
determinants  of  the  appropriate  principal  minors,  as  predicted  by  (14)  and  (15). 

Although  there  has  been  some  growth  in  the  magnitudes  of  the  matrbc  elements  it 
has  been  kept  in  much  better  check.  The  whole  of  this  computation  could  have  been 
achieved  using  standard  16-bit  integer  arithmetic  -  even  including  the  temporary  values. 
The  division-free  algorithm  in  this  case  require  at  least  20  bits  even  though  the  original 
matrix  has  integer  elements  bounded  by  8. 

Remark  17.  It  is  necessary  to  observe  that  the  vaiious  fraction-free  versions  of  GE  are 
well-suited  to  conventional  binary  integer  arithmetic  -but  cannot  easily  be  applied  hi  other 
integer  arithmetic  implementations  such  as  RNS  since  the  divisions  can  not  be  peiformed 
within  the  RNS  system  itself.  This  is  true  even  in  the  ideal  situation  of  Algorithm  14 
where  all  divisions  are  known  to  have  exact  integer  results. 
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3.4.  Effect  on  computation.  In  this  section,  we  describe  the  effect 
fraction-free  GE  Algorithm  14  on  the  solution  of  the  various  underlying  pr  • 

taL  U«  action,  it  is  Plata  that  the  ^aation  .s  par.K=ul.rly 

simple,  provided  only  that  no  zero  pivots  arise  during  the  computatio  . 

Linear  systems.  For  the  solution  of  a  (nonsingular)  linear  system  using  Algorithm 
14  there  is  o^  course  no  guarantee  that  the  solution  vector  has  only  integer  components. 
lUhe^sL  is  known  (b^nsa  of  the  context,  perhaps)  to  have  an  mteser  solutjon  then 
this  can  be  computed  by  applying  the  conventional  back  substitution  (Algorithm 
Algorithm  8,  for  Sample)  in  which  case  the  divisions  will  again  have 
Otherwise  the  solution  can  be  express«i  as  fractions  by 
In  its  column  form  equivalent  to  Algorithm  8,  the  rational  arithmetic  back  sul^itution 
can  be  described  as  follows.  The  algorithm  below  makes  use  of  a  function  reduce  [p,  9] 
which  simply  reduces  a  rational  number  p/g  to  its  normal  form  by  removing  any  common 

factors. 

Algorithm  15  Back  substitution  using  rational  arithmetic  —  column  form 
Input  n  X  n  upper  triangular  matrbc  U,  n-vector  ^ 

IniMaUze  mtional  solution  m  :=  0;  n  :=  1  (solution  ml\  be  =  m  rii) 

Rationalize  right-hand  side  p  :=  b;  q  :=  1  (right-hand  side  6,  =  Pi/gi) 

Compute 

for  j  =  n  downto  1 
rrij  :=Pf  rij  := 
reduce  [mj ,  71  j] 
for  r  =  1  to  j  —  1 

Pi  pi*  77. j  —  Uij  *  777. j  *  Qi 

Qi  :=  qi  *  nj 
reduce  [pi,  qi] 

Output  solution 

Remark  18.  Note  that  the  inner  loop  call  to  reduce  [p,  ,9,]  could  be  eliminated  at  little 
expense  It’s  removal  would  have  the  benefit  that  the  same  denominator  would  be  used 
for  each  bi  for  i  <  j.  The  cost  in  terms  of  dynamic  range  would  likely  be  mmimal  since 
Z%ted  range  would  only  be  decreased  if  everj-  call  to  reduce b. : some  value 
of  j  results  in  a  smaller  value  of  qi. 

Example  4.  Consider  the  coefficient  matrix  of  Example  3  with  the  original  right-hand 
side  vector  [25,20, 25,20]^. 

Algorithm  14  would  generate  the  augmented  matrbc 


■  8 

7 

4 

1 

25 

0 

20 

40 

20 

60 

0 

0 

no 

150 

260 

0 

0 

0 

—450 

—450 

Applying  Algorithm  15.  we  first  get,  for  j  =  4.  7714  =  -450, 774  =  -|50  ®»^^® 
of  reduce  [7774,774]  is  then  n7.4  =  774  =  L  The  inner  loop,  with  the  toher  c^l  o_iedu^ 
then  yields  p,  =  24, P2  =  40.p3  =  HO  and  94  =  92  =  93  =  I-  Next-  "ith  .j  -  -3. 
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get  m3  =  n3  =  1  and  then  pi  =  20,P2  =  0  with,  still,  9i  =  92  =  1-  The  next  iteration 
of  the  i-loop  gives  m2  =  0,n2  =  l,Pi  =  20  and  91  =  1  from  which,  finally,  we  get 
rrii  =  20,  ni  =  8  which  reduces  to  mi  =  5,  n-i  =  2  so  that  the  complete  rational  soluti^  is 
[5/2, 0/1, 1/1,  l/l]^-  In  this  particular  case  the  reduction  in  the  inner  loop,  had  no  effect 
since  the  outer  loop  reduction  resulted  in  ^2  =  ^  =  94  =  I* 

Determinant.  We  have  already  seen  that  Algorithm  14  delivers  det  A  automatically 
as  the  final  value  of 

Rank  determination.  In  the  absence  of  a  zero  pivot,  again  there  is  nothing  further 
to  be  done.  In  the  event  of  such  a  zero,  then  interchanging  the  pivot  row  with  any  lower 
row  with  a  nonzero  entry  in  the  pivot  column  will  allow  the  fraction-free  algorithm  to 
proceed.  Simply  counting  the  number  of  nonzero  entries  on  the  diagonal  gives  the  rank, 
€is  before.  This  simple  interchange  is  not  necessarily  the  optimal  pivoting  strategy  for 
integer  computation. 

In  the  event  of  a  rank-deficient  matrix,  obtaining  the  solution  space  for  an  underde¬ 
termined  system  can  be  accomplished  by  modifying  Algorithm  15  in  just  the  same  way 
as  was  proposed  in  Section  2.3  for  real-number  computation. 

3.5.  Pivoting.  Clearly  in  the  event  of  a  zero  pivot  some  pivoting  is  necessary  in 
any  good  implementation  of  GE.  The  question  is  w’hat  is  the  right  strategy  for  integer 
computing?  Choosing  the  largest  element  in  the  floating-point  algorithm  has  the  virtue  of 
keeping  all  multipliers  small  and  therefore  restricting  the  growth  of  the  matrix  elements. 
However  that  restricted  growth  is  only  realized  because  of  the  divisions  that  are  performed 
in  Algorithms  1  through  6. 

At  least  intuitively,  choosing  the  largest  element  in  the  pivot  column  is  not  necessarily 
good  for  integer  computation.  Indeed  a  large  prime  pivot  is  probably  near-worst  since 
that  appears  to  almost  guarantee  rapid  growth  in  the  dynamic  range.  In  (11)  and  then 
(12)  we  see  that  each  has  the  factor  an  and  each  further  stage  has  this  factor  and 

then  has  it  repeated  in  the  factors  The-earlier  a  particular  pivot  is  used  the  greater 

the  impact  it  has  on  subsequent  range  gro\nh.  This  is  made  plain  in  equation  (13). 

In  the  fraction-free  version  Algorithm  14.  however,  the  range  gro\^iih  is  essentially 
independent  of  the  order  of  the  use  of  the  rows.  In  the  case  of  a  full  rank  matrix  A, 
the  final  entry  qnn  —  ^  (I^)  ^  necessarily  the  largest  of  the  principal  minor 

determinants  and  this  v-alue  is  invariant  (up  to  sign  changes)  under  row  interchanges.  It 
follows  that  the  simplest  pivoting  strateg>’  is  probably  the  best  for  this  algorithm.  That 
is,  at  stage  we  should  use  the  first  row  for  which  the  pivot  column  entr}^  is  nonzero: 
choose  pivot  =  min  {p  >  ?'  :  Up,  0}. 

It  is  concei\^ble  that  a  different  ordering  may  represent  some  minor  improvement  on 
this.  For  example,  looking  for  pivots  which  are  either  powers  of  the  binary  (or  other  com¬ 
putational)  base  may  make  subsequent  divisions  particularly  simple.  However,  searching 
for  the  appropriate  pivot  element  in  this  regard  would  (almost  surely)  be  moie  wasteful 
than  beneficial. 

3.6.  Complexity  revisited.  In  this  section,  we  begin  by  obtaining  the  basic  integer 
operation  counts  for  the  fraction-free  GE  Algorithm  14.  There  are  two  essential  differences 
between  these  counts  and  those  for  the  floating-point  algorithms  presented  in  Section 
2.4:  the  fundamental  elimination  loops  contain  no  divisions  but  ha\e  twice  as  man^> 
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multiplications,  and  there  is  the  additional  complexity  of  the  removal  of  the  common 

^^'^he^additkmarmultipliciions  in  the  elimination  are  easily  counted.  ^ 
factor  removal  entails  only  divisions.  The  total  number  of  these  can  be  assessed  from  the 
For  eL.  i  >  2,  ther.  U  one  division 

are  solvine  a  system)  which  runs  from  j  =  t  + 1  to  n.  Also  in  this  loop  is  the  l-top  w  hich 
has  the  same  limits  knd  contains  another  division.  The  total  division  connt  is  therefor. 

-  i)(n  -  r  +  1)  =  53 !(/  +  1)  =  -^^{71  -  l)(n  -  2) 


The  total  operation  counts  are  summarized  in  Table  4. 

TABLE  4  Integer  arithmetic  operation  counts  for  the  fraction-free  GE  Algorithm  14  for 
solution  of  an  n  x  u  linear  system  -4x  —  b. 


Operation  Matrbc  Factorization 


-h  or  - 


burr  -  1 


|n(n2-l)  n(n-l) 

l(n  -  2)  (n  -  1)  (2n  -  3)  ^(n  -  l)(n  -2) 


Right*hand  side  TOTAL _ 

^n{n-lY~  |n(n-l)(2n  +  5) 


|n  {n  —  1)  (2n  +  5) 
|n  (n  —  1)  (n  —  2) 


Note  that  this  operation  count  does  not  include  the  back  substitution  Algontom  lo. 
What  we  see  most  importantly  is  that  the  number  of  multiplications  has  donh^  a^^ 
even  worse,  the  number  of  divisions  has  increased  by  a  complete  order  from  0(n  )  to 

^^Unfortunately  this  is  not  the  end  of  the  story.  The  complexity  of  this  algorithm  is 
further  increased  by  virtue  of  the  fact  that  these  integer  divisions  are  typically  inore 
difficult  than  their  floating-point  counterparts  -  especially  with 

we  have  already  seen  can  be  substantial.  This  is  likely  to  necessitate  the  use  of  long 
integer  arithmetic  for  which  special  algorithms  must  be  used. 

3  7.  Arithmetic  with  long  integers.  In  this  section  we  discuss  briefly  some  aspects 
of  the  long  integer  arithmetic  which  is  required  for  integer  Gauss  elimination  alS°"thms. 
Long  integer  arithmetic  can  be  simulated  using  multiple  words  ^  °  ' 

For  exlmple.  if  the  underlying  integer  wordlength  is  8  bits  (or  1  bvte)  then  such  an 
integer  can  be  regarded  as  a  radix  2«  =  256  digit.  Conventionally  the  range  of  ^®lues  of  the 
base^256  digits  would  be  -128,-127, . . . ,  127  so  that  signed  integers  can  be  represented. 
For  simplicity  in  the  current  description,  restrict  our  attention  to  nonnegative  integers 
with  a  digit  range  of  0, 1, ,  255.  Very  large  integers  could  then  be  stored  using  a  vector 

of  such  integers  using  conventional  place  value.  _  , 

In  general,  suppose  the  base  wordlength  is  I  bits  so  that  the  effective  rad^  «  R- 2 

The  vector  (d<,,di,...  ,dK-i)  would  then  be  used  to  represent  the  integer  A  €  [O.R  IJ 

y  =  •¥dK-2R^~^ - l-difi-fdo 

where  each  digit  satisfies  0  <  d,  <  R.  For  efficient  arithmetic  using  such  a  representa¬ 
tion  we  require  an  integer  accumulator  with  2L  bits.  Among  other  consequence  of  this 
are  that  Sdition  can  be  computed  in  “digit-paraller’  with  subsequent  attention  to  am^ 
carrie  which  mav  be  propagated.  A  large  radix  carry  lookahead  or  conditional  sum  adder 
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could  be  constructed  to  improve  the  efficiency  of  addition.  These  addition  algorithms  are 
simple  generalizations  of  the  usual  binary  addition  algorithms  which  are  described,  for 

example,  in  [19].  .  u 

Alternative  approaches  to  long  wordlength  integer  arithmetic  are  provided  by  using 
Residue  Number  Systems,  RNS.  This  is  a  naturally  parallel  integer  arithmetic  which  does 
not  readily  admit  division.  It  is  therefore  not  very  suitable  for  use  with  the  fraction-free 
Gauss  elimination.  Algorithm  14.  It  has  been  considered  in  detail  for  naval  applications 
in  [9],  [10]  and  [11]. 

Convolution  form  of  product.  Multiplication  of  long  integers  can  be  achieved 
with  reasonable  efficiency  using  the  convolution  form  of  the  product  working  with  the 
component  words  of  the  large  integers.  Again  for  simplicity,  we  shall  only  consider  mul- 
tiplication  of  positive  integers.  The  multiplication  of  two  integers  in  the  form  (16)  will 
result  in  an  integer  whose  representation  requires  at  most  2K  L-bit  words. 

Suppose  that  we  require  the  product  of  the  long  integers 

K-l 

7n=^Qii?\  n  = 
issO 


This  product  is  given  by 

/K-i  \  /a'-i  \  2A'-2/^l  \ 

(17) 

\t=o  /  \;-0  J  \  t-o  / 


where  we  use  the  convention  I3j  =  0  if  j  <  0.  This  is  essentially  a  convolution  product  of 
the  coefficient  vectors. 

Now  each  coefficient  product  aipk^i  can  be  viewed  as  a  hase-R  digit  which  we  can 
write  in  the  form  aij3k-i  =  +  where  each  7/e,t  and  6k^i  is  a  hase-R  digit.  Thus 

is  the  most  significant  and  6k^i  the  least  significant  i?-digit  of  ailSu^i.  The  product 
(17)  can  therefore  be  written  as 


2K-2  /K-l 


m  * 


fc=0  \  t:=0 


2K-1  /K-l 


(18) 


where  7^,1  =  ^k,i  =  0  for  i  >  k. 

Remark  19.  Carries  beyond  the  position  is  not  possible  since  m*n  <  (R^-1)^  < 

Remark  20.  The  results  of  the  inner  sums  in  (18)  will  usually  create  further  carries 
which  must  be  accounted  for.  However,  we  note  that  for  almost  any  minima]  degiee  of 
parallelism  in  the  processor,  each  of  these  inner  sums  can  be  performed  simultaneously 
since  they  consist  of  at  most  2J\  terms  each  less  than  R  and  we  would  expect  2I\  R. 
Therefore  the  sum  will  be  less  than  R^  so  that  it  will  be  representable  in  a  double  length 

accumulator. 

For  the  8-bit  basic  wordlength  the  restriction  is  only  that  our  integers  do  not  exceed 
>  10®^^  which  is  well  beyond  any  topical  integer  computing  range  for 
linear  algebra  problems. 
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The  length  of  the  inner  sums  in  (18)  will  also  restrict  the  size  of  any  carry  and  therefore 
may  be  useful  in  bounding  the  range  of  the  propagation  of  such  carries.  This  ^uld  be 
used  to  improve  the  efficiency  of  such  a  convolution  product.  We  do  not  consider  such 

details  further  in  this  report.  i  •* 

What  effect  does  such  a  long  multiplication  have  on  the  arithmetic  complexi  y 
an  integer  algorithm?  The  product  formula  (17)  entails  K  basic  multiplications,  the^ 
components  must  then  be  broken  into  their  two  component  digits  and  then  approximately 
A*  additions  together  with  the  carries  these  generate.  In  Table  4,  th^  means  that  each 
of  the  (approximately)  |n®  multiplications  entails  something  like  2A  regular  integer 
arithmetic  operations.  For  even  quite  moderate  range  growth  with  A  -  4  this  has  the 
effect  of  increasing  this  part  of  the  complexity  by  a  factor  of  32. 


Division.  Division  of  long  integers  can  also  be  achieved  by  generalizing  some  of  the 
standard  algorithms  which  are  used  in  biiiarj-  integer  hardware  such  as  the  SRT  division 
algorithm  which  is  based  on  the  idea  of  nonrestoring  division.  This  algorithm  is  well-suited 
to  high-radix  division  and  so  can  be  modified  to  the  long  integer  framework. 

The  SRT  algorithm  relies  on  repeated  addition  and  subtraction  using  signed  digits. 
Since  we  have  only  dealt  with  arithmetic  of  nonnegative  long  integers  here,  we  do  not 
discuss  the  detailed  implementation  further.  For  details  of  the  basic  algorithm  and  its 
implementation  at  least  for  radix  4,  see  [19]. 


3.8.  Should  we  use  real  arithmetic  anyway?.  Unless  roundoff  errors  are  severe, 
it  may  well  be  prudent  to  compute  the  solutions  to  integer  linear  algebra  problems  b\ 
applying  the  appropriate  real,  floating-point,  algorithm  and  rounding  the  results.  All  ^e 
difficulties  of  range  growth,  avoiding  divisions  or  additional  complexity  are  then  removed. 

Ftom  the  error  analysis  summarized  in  Section  2.5,  we  can  obtain  estimates  of  the 
number  of  correct  digits  in  the  final  result.  First  note  that  since  the  underlying  pioblem 
has  an  all-integer  coefficient  matrix  (and  right-hand  side)  then  we  may  assume  that  we 
have  exact  data.  The  only  contribution  to  the  error  in  the  solution  of  such  a  linear  system 
therefore  arises  from  the’  roundoff  errors  of  floating-point  arithmetic.  Only  foi  very  ill- 
conditioned  matrices  is  the  upper  bound  on  the  error  given  by  the  right-hand  side  of  the 
inequality  (9)  likely  to  indicate  that  the  rounded  solution  of  an  integer  system  may  be 
incorrect.  With  the  improved  elementwise  bounds  obtained  by  Demmel  and  others  even 
greater  confidence  in  the  computed  results  would  be  obtained.  (See  [4]  for  an  introduction 
to  the  literature  on  this  subject.)  .  _  .  .  j  w 

We  have  also  seen  that  the  complexity  of  integer-arithmetic  GE  is  increased  relative 
to  the  floating-point  algorithms.  Compare  Tables  2  and  4.  By  its  nature  the  fraction- 
free  Algorithm  14  is  “row-oriented”  which  may  render  it  advantageous  on  some  parallel 
architectures,  for  instance  the  removal  of  the  common  factors  would  be  a  simple  array- 
operation  on  a  massively  parallel  array  processor.  However  whatever  advantages  could  be 
realized  for  this  algorithm  could  also  be  obtained  for  the  appropriate  arrangement  of  the 
original  LU  factorization  algorithm  such  as  Algorithms  4,  5  or  6. 

On  most  modern  computer  architectures  the  speed  of  floating-point  operations  is  com¬ 
parable  with  all  but  the  shortest  of  integer  wordlengths.  The  range  growth  problem  of 
integer  computation  makes  such  short  wordlengths  impractical  for  most  problems^  If, 
using  floating-point,  the  arithmetic  is  essentially  no  slower,  the  complexity  is  reduced  be¬ 
cause  division  is  available  and  the  accuracy  of  integer  arithmetic  can  be  recovered  in  the 
final  results,  it  seems  this  (real  arithmetic)  would  be  better  in  most  circumstances. 
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4.  Over  the  ration als 

In  this  section, we  consider  the  use  of  Gauss  elimination  in  the  setting  of  rational  arith¬ 
metic  Of  course  in  some  sense,  floating-point  is  rational  arithmetic  but  it  uses  a  ver>' 
special  subset  of  the  rationals  and  does  not  comply  with  the  axioms  of  conventional  ratio¬ 
nal  arithmetic.  We  are  concerned  here  with  matrices  with  entries  in  the  field  Q  of  rational 
numbers.  The  arithmetic  operations  will  be  similar  in  nature  to  those  of  Algorithm  15  for 
back  substitution  in  the  integer  GE  solution  of  a  linear  system. 


4.1.  Rational  representations.  There  are  choices  to  be  made  over  the  way  in  which 
rational  numbers  are  to  be  stored  and  manipulated  within  the  computer.  The  conceptually 
simplest  option  is  simply  to  store  9  €  Q  as  a  pair  of  integers  representing  its  numerator 
and  (positive)  denominator  in  its  maximally  reduced  form.  Alternatives  that  have  been 
extensively  resecurched  include  the  use  of  continued  fraction  representations  and,  in  par¬ 
ticular,  the  lexicographic  continued  fraction,  or  LCF,  arithmetic  of  Kornerup  and  Matula 
[13],  [14],  [15],  for  example. 

[n/m]  form.  The  standard  representation  of  rational  numbers  is  as  a  quotient  of  two 
co-prime  integers  ^ 

^  "  m 

where  n,m  €  Z  have  no  common  factors  and  m  >  0.  Arithmetic  operations  are  then 
defined  according  to  the  usual  rules  of  rational  arithmetic  with  reduction  to  this  “normal¬ 
ized”  form  after  each  arithmetic  operation. 

It  is  immediately  apparent  that  many  of  the  range  g^o^^^h  problems  which  plague 
integer  GE  will  reapp>ear  here  with  comparable  severity.  Of  course  the  divisions  w’hich  are 
inherent  in  the  basic  forms  of  GE  and  LU  factorization  can  be  performed  here  and  restrict 
the  range  of  values  of  the  rational  numbers  being  represented  -  but  these  divisions  do  not 
necessarily  reduce  the  range  of  integers  needed  for  the  numerators  and  denominators 
sep)arately. 

The  chief  virtue  that  would  be  derived  here  is  the  elimination  of  any  rounding  errors 
and  therefore  definitive  answers  to  questions  such  as  the  rank  of  the  matrix  and  exact 
values  for  the  determinant  and  for  the  solution  vector  of  a  linear  system..  This  last  would 
be  computed  using  the  (obvious  and  minor)  modification  of  Algorithm  15. 

Use  of  continued  fractions.  An  alternative  representation  of  rational  numbers  is 
provided  by  using  continued  fractions.  This  possibility  has  been  extensively  explored  by 
Kornerup  and  Matula  in  a  substantial  series  of  papers  including  [13],  [14]  and  [15].  In  the 
later  papers  in  this  series,  techniques  are  described  for  performing  rational  arithmetic  in 
a  bit-pipelined  manner.  This  has  the  efTect  of  allowing  many  simple  rational  functions  of 
arguments  expressed  in  continued  fraction  form  to  be  computed  much  more  quickly  than 
might  otherwise  be  the  case.  This  is  based  in  p>art  on  the  use  of  a  particular  encoding  of 
the  numbers  which  has  retains  the  natural  ordering  of  the  representations.  This  is  called 
the  lexicographic  continued  fraction  representation  or  LCF.  Its  details  are  not  important 

here.  ^  . 

The  continued  fraction  representation  of  a  positive  rational  number  is  given  by 


g  =  00  + 


1 


where  the  integers  oo  >  0,  and  o,  >  1  for  ?  >  1- 


1 


(19) 
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The  manifesution  of  the  range  growth  problem  here  hes  not  m  the 
the  various  coefficients  so  much  as  in  the  mmber  of  them;  th^  of  LCF 

strine  oo  a,  02....  used  to  represent  a  particular  rational  num^.  For  details  of  LCF 
arithmSc,  its  implementation  and  its  use  the  reader  is  referred  to  the  original  pape^ 
This  section  is  included  solely  to  increase  awareness  of  the  existence  of  alternative  for 
of  rational  arithmetic  which  could  be  used. 

4.2.  GE  algorithm  for  rational  arithmetic.  The  algorithm  for  Perming  the 
factorization  of  a  rational  matrix  can  be  any  of  the  variants  disci^  in  Secti^  2  «  »th 
the  real  arithmetic  operations  replaced  by  rational  arithmetic^  There  is  no 
deteUing  all  these  algorithms  explicitly.  In  order  to  disc^  the  dynamic  rang  ^ 
this  context  however  it  is  helpful  to  have  at  least  the  basic  ijfc-form  described  for  this 
context.  We  shall  denote  the  original  matrix  elements  by 

P«i 

and  those  of  the  right-hand  side  vector  by 

where  all  fractions  are  assumed  to  be  in  reduced  form. 

Algorithm  16  Basic  GE  for  rational  arithmetic  ijk  form 

Input  nxn  matrix  rational  A  (and  right-hand  side  b  if  solving  a  system) 

Compute 

for  i  =  1  to  n  —  1 
for  j  i  +  1  to  n 

TTlji  :=  ^ji  •=  Pii^ji 

reduce  [mj<,nji) 

Tj  :=  rjUjiSi  -  Sjmjiri;  s,  =  sjUjiSi  (for  system) 
reduce  [rj,Sj] 
for  fc  =  i  -I- 1  to  n 

pjk  :=  Pik^ji^ik  -  qjkmjiPik',  <ljk  =  gjknjiqik 
reduce  (pjjr.flji.] 

Pji  :=  0;  Qji  .—  1 

Output  (modified)  matrbc  A  (and  b  if  solving  a  system) 

Remark  21.  Algorithm  16  is  just  the  basic  Algorithm  1  modified  for  rational  arithmetic. 
It  is  plain  that  the  arithmetic  complexity  has  been  increased  considerably  as  a  result  of 
storing  each  element  as  a  quotient  of  integers  and  by  the  need  for  reduction  of  these 
rational  numbeis  to  their  “normalized''  form. 

No  pivoting  has  been  incorporated  here  and  as  with  the  standard  integer  setting  it  is 
not  obvious  what  the  optimal  pivoting  strateg>-  would  be.  In  the  above  algorithm  it  is 
simplv  assumed  that  no  zero  pivots  are  encountered.  If  such  a  pivot  element  does  in  fac 
arise  we  can  simply  interchange  the  piwt  row  with  a  row  which  has  a  nonzero  entr>  in 

the  pivot  column. 
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From  the  point  of  view  of  range  growth  there  may  be  some  merit  in  choosing  pivote  to 
be  as  “simple”  as  possible.  Simple  here  would  correspond  roughly  to  haying  the  small«t 
denominator  (or  maybe  the  denominator  whose  largest  prime  factor  is  minimum).  Such  a 
search  for  the  simplest  pivot  may  well  be  more  expensive  than  is  merited  by  ariy  reduction 
which  may  be  achieved  in  the  dynamic  range.  The  conventional  argument  for  pivoting 
(apart  from  avoiding  zero  pivots)  is  not  relevant  to  either  integer  or  rational  arithmetic 
since  roundoff  error,  and  the  associated  question  of  numerical  stability,  are  not  concerns. 

4.3.  Complexity.  The  most  obvious  increase  in  the  complexity  in  Algorithm  16  arises 
out  of  the  mere  fact  that  it  is  performing  rational  arithmetic  and  so  every  arithmetic 
operation  entails  manipulation  of  both  the  numerator  and  denominator.  There  is  also 
the  further  complication  arising  from  the  reduction  of  each  ordered  pau-  to  represent 
an  irreducible  fraction.  This  requires  either  that  each  integer  is  stored  as  a  product 
of  its  prime  factors,  or  more  reasonably,  that  the  Euclidean  algorithm,  Algorithiri  12, 
for  frnding  the  gcd  of  two  integers  is  employed  and  followed  with  two  integer  divisions. 
Because  the  Euclidean  algorithm  is  iterative  we  cannot  determine  the  number  of  integer 
mod  operations  that  are  required.  (Also  any  bounds  which  could  be  derived  would  be 

hopelessly  pessimistic.)  *  ^ 

In  Table  5,  we  list  the  numbers  of  conventional  integer  arithmetic  operations  together 
with  the  number  of  gcd’s  that  are  needed.  Typically  we  might  expect  that  a  gcd  would 
be  equivalent  to  several  integer  divisions  and  that  these  divisions  are,  in  turn,  equivalent 
to  several  multiplications.  However  it  should  also  be  noted  that,  because  of  the  range 
growth,  the  multiplications  may  involve  long  wordlength  integers.  On  the  other  hand 
we  may  expect  the  gcd  to  be  much  smaller  so  that  the  divisor  wordlengths  may  be  less 
extreme.  In  the  operation  counts  in  Table  5,  no  attempt  is  made  to  account  for  dynamic 
range  considerations  or  the  relative  weights  to  be  given  to  the  different  operations. 


TABLE  5  Integer  arithmetic  operation  counts  for  the  rational  GE  Algorithm  16  for 
solution  of  an  n  X  n  linear  system  Ax  =  b. 


Operation 

Matrix  Factorization 

Right-hand  side 

TOTAL 

+  or  - 

X 

/ 

gcd 

(n^  -  1) 

2n^  (n  —  1) 

in(n-  1) 

3n(n  —  1) 
n(n  —  1) 
in{n-  1) 

in(n  -  1)  (2n  +  5) 
n  (ti  —  1)  (2n  +  3) 
in(n  —  l)(2n-i-5) 
|n  (n  —  1)  (2n  +  5) 

Again  the  biggest  single  effect  is  that  the  number  of  divisions  has  increased  to  0(n  ) 
in  addition  to  the  O  [n^)  gcd  operations.  This  time  the  number  of  multiplications  has 
also  increased  by  a  substantial  factor.  The  overall  operation  counts  suggest  that  the  cost 
of  rational  computation  may  be  too  high  to  justify  the  increased  stability  and  accuracy 
especially  when  the  range  growth  and  resulting  additional  complexity  of  these  operations 
is  considered  too. 

4.4.  Growth  in  dynamic  range.  In  attempting  to  analyse  the  potential  growth  in 
the  necessary  dynamic  range  for  implementing  Algorithm  16,  we  cannot  assume  any  useful 
reduction  in  the  rational  quantities  other  than  that  which  is  a  necessary  consequence  of  the 
elimination  procedure  itself.  For  simplicity  in  this  setting,  we  make  the  assumption  that 
the  original  matrix  A  is  in  fact  an  integer  matrix.  This  allows  us  to  compare  the  results 
of  the  (rational)  Algorithm  16  with  those  that  would  be  obtained  using  the  divisionless 
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integer  version  Algorithm  11.  In  much  the  same  way  as  was  described  for  the  fracton-free 
Algorithm  14,  we  shall  see  that  this  comparison  reveals  certain  comnmn  factors  '^hich  v  il 
be  removed  by  the  various  reduction  steps.  This  has  the  effect  of  -educing  the  potential 
range  growth  in  a  similar  manner  to  that  which  we  observed  in  Algorithm  14. 

As  for  the  analvsis  in  Section  3.3,  we  shall  concentrate  on  the  matrix  iteelf  Similar 
analysis  is  ^’alid  for  the  right-hand  side  of  a  linear  system  but  ^ds  nothing  to  overall 
picture  Again  following  the  earlier  analysis,  we  shall  denote  by  Oy  the  mteger  elements 
of  the  origiMl  matrix.  The  results  of  the  i-th  stage  of  the  integer  Algorrthm  ll^will^in 
be  denoted  by  The  corresponding  rational  values  will  be  denoted  by  Cj^.  These 
are  of  course  just  rational  representations  of  the  quantities  generated  by  the  original 

the  first  stage  of  the  elimination,  the  multipliers  used  are  (irreducible  rational 
representations)  of  ^ 

TTlji  = 

ail 

Except  in  the  unlikely  event  that  all  elements  in  the  first  column  of  A  share  a  wmmon 
factor,  the  range  for  the  denominator  cannot  be  reduced.  The  results  of  the  first  stage  are 

(1)  fljl  QiiQflr-  QjlOlfc  _  (20) 

on  ■  ‘111 

and  again  unless  there  is  some  general  cancellation  the  range  growth  for  the  numerator  is 
identical  to  that  for  Algorithm  11  at  this  stage. 

Similarly  the  second  stage  results  in 

M) 

J2)  _  Jl)  _  filciV 
^22  , 

and  we  observe  that  in  light  of  (20)  the  multiplier  here  is  the  same  as  would  be  used  with 
the  integer  matrix.  Indeed,  substituting  (20)  into  the  last  equation  we  get 


.(2)  _  ^  _  ^l2^fc  =  °22  <’ik  -  Pj2  '^2k  ^ 


■  a„  60)au 


^22^  ttl  1 


b^  Ol  1 


However,  we  recall  from  Section  3.3  that  each  bf^  has  the  factor  an  which  will  Jerefore 
be  removed  by  the  rational  reduction  step.  That  is  the  only  guaranteed  factor  which  can 
be  used  to  help  control  the  dynamic  range  growth.  The  range  required  for  teh  numerator 
at  this  stage  is  therefore  precisely  the  same  as  that  for  teh  fraction-free  Algorithm  1  . 
The  denominator  (assuming  no  further  general  common  factor  is  present)  is  622  -  <^2,  the 

determinant  of  the  2x2  principal  minor.  •  i  + 

Because  of  this  common  denominator,  the  next  stage  is  again  essentially  equivalent 

to  that  for  Algorithm  14  since  becomes  a  common  factor  in  the  remaining  actn'e 
matrix.  Ultimately  then  the  range  growth  for  Algorithm  16  is  identical  to  that  which 
is  encountered  in  the  fraction-free  Algorithm  14.  The  inference  to  be  drawn  from  this 
analysis  is  that  if  the  initial  matrix  is  an  integer  matrix,  using  rational  arithmetic  jields 
no  benefit  relative  to  the  fraction-free  algorithm.  ,  .  , 

If  the  original  matrix  consists  of  rational  entries,  the  range  growth  problem  become 
potentially  even  more  critical  since,  the  “cross-multiplications’  needed  for  the  elimination 
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do  not  appear  to  generate  any  obvious  and  general  common  factors.  Unless  the  original 
denominators  are  small  and  have  similar  factors  it  seems  unlikely  that  any  purely  rational 
algorithm  will  be  practical  and  the  best  solution  would  almost  certainly  be  to  use  real 
(floating-point)  arithmetic. 


5.  General  rings 

In  this  section,  we  consider  the  applicability  (and  application)  of  GE  in  a  more  general 
algebraic  setting.  There  are  at  least  two  fundamental  questions  which  arise  immediately: 

“When  do  the  problems  have  solutions?’  and 

“When  do  the  algorithms  make  sense?” 

We  are  interested  here  in  matrices  with  entries  in  a  general  ring  Tl.  In  order  to  answer 
the  fundamental  questions,  we  shall  need  to  impose  further  conditions  on  this  ring.  For 
details  of  the  definitions  and  properties  of  the  various  algebraic  structures  see  [5]  or  any 
standard  text  on  abstract  algebra.  The  key  properties  which  are  being  used  will  be 
summarized  as  they  are  introduced  here. 

In  essence  a  ring  72.  is  a  structure  which  has  both  “addition”  and  “multiplication” 
defined  on  it.  The  addition  has  all  the  properties  associated  with  the  various  number 
systems.  The  multiplication  need  not  have  all  the  properties  of  the  real  number  system. 
In  a  general  ring,  all  we  may  assume  is  that  the  multiplication  is  associative  and  that  it 
is  distributive  over  addition.  “Vector  spaces”  over  rings  are  called  modules.  The  concepts 
of  linear  independence,  basis  and  scalar  multiplication  carry  over  in  a  natural  way  to 
modules. 

Special  cases  of  rings  include  fields  of  which  the  rational,  real  and  complex  numbers 
Q,  R  and  C  are  the  most  important  examples,  and  integral  domains  such  as  the  integers 
Z  and  the  ring  of  polynomials  over  the  real  numbers  R[.t].  In  a  field,  multiplication  and 
division  (except  by  zero)  are  defined  and  have  the  expected  properties.  In  an  integral 
domain,  the  multiplication  is  commutative,  there  is  a  multiplicative  identity,  1,  and  there 
are  no  divisors  of  zero.  The  integers  themselves  have  further  properties  which  are  relevant 
to  the  above  questions.  Specifically,  the  integers  form  a  unique  factorization  domain 
which  is  to  say  that  every  positive  integer  n  7^  0, 1  can  be  written  as  the  product  of  prime 
numbers. 

Associated  with  every  integral  domain  is  its  field  of  fractions.  For  the  integers  this 
field  is  just  Q,  the  rational  numbers.  The  field  of  fractions  Q  for  a  general  integral  domain 
.Z  is  a  field  which  results  from  the  analogous  construction  and  consists  of  ordered  pairs 
(m,n)  where  m,n  €  Z  with  its  “arithmetic”  defined  just  as  for  real  fractions.  As  was  the 
case  in  Section  3.4,  it  may  be  the  case  that  a  linear  system  of  equations  has  no  solution 
in  2  but  does  have  a  “rational”  solution  in  its  field  of  fractions  Q. 

5.1.  Problem  definitions.  Before  discussing  any  algorithms  or  the  algebraic  struc¬ 
tures  within  which  they  make  sense,  it  is  necessary  to  reconsider  the  definitions  of  the 
various  linear  problems  under  discussion  and  to  determine  when  the  questions  themselves 
are  well-defined. 

Since  the  concept  of  linear  dependence  carries  over  to  arbitrary  modules,  we  may 
define  the  rank  of  a  matrix  A  with  entries  in  a  general  ring  Tl  by  the  maximum  number 
of  linearly  independent  rows  of  A.  The  question  of  determining  the  rank  of  A  therefore 
can  be  addressed  for  a  completely  general  ring  72.  We  shall  also  see  that  using  (a  minor 
modification  of  the  divisionless)  Algorithm  11  will  solve  this  problem. 
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Simil&rlv  the  determinant  can  be  defined  for  matrices  over  a  general  ring  using  the 
“permutation*based”  definition  for  the  expansion  of  a  determinant.  However,  the  divi¬ 
sionless  GE  algorithm  does  not  readily  yield  the  detA  even  for  the  integers  Z  due  to  the 
growth  of  the  elements  resulting  from  the  elementary  row  operations  performed.  More 
structure  is  needed  in  the  ring  7^  if  we  are  to  be  able  to  obtain  det-4  using  GE.  The 
fraction-free  Algorithm  14  can  be  performed  in  any  unique  factorization  domain.  The  ba¬ 
sic  GE  Algorithm  1  can  be  modified  to  solve  the  problem  in  a  (noncommutative)  division 
ring. 

Something  similar  to  division  is  clearly  needed  in  order  that  a  system  of  linear  equations 
can  be  solved.  Even  when  such  a  system  is  well-defined,  it  may  be  necessary  to  go  to  the 
field  of  fractions  in  order  to  obtain  the  solution.  So  for  the  solution  of  systems  of  linear 
equations,  the  underlying  algebraic  ring  must  have  most  of  the  properties  of  the  integers 
and  the  familiar  number  systems. 

5.2*  Algorithms. 

Rank.  The  division-free  GE  Algorithm  11  is  well-defined  for  any  commutative  ring. 
For  a  completely  general  ring  the  only  modification  needed  is  that  the  order  of  the 
various  products  must  be  consistent.  The  notion  of  rank  is  also  well-defined.  The  following 
algorithm  is  a  simple  modification  of  Algorithm  7. 

Algorithm  17  Division- free  GE  in  a  general  ring  Tl 

Input  n  xn  matrix  A  with  entries  from  72, 

Compute 

for  7  =  1  to  n  —  1 
for  j  =  z  -h  1  to  n 
for  fc  =  t  -j- 1  to  71 

Q>jk  Q-U-Q'jt 

Qji  :=  0 

Output  (modified)  matrix  A 

Remark  22.  The  order  of  the  products  in  the  elimination  step  is  chosen  to  be  consistent 
j\’ith  the  elimination  of  aji. 

Remark  23.  In  order  to  obtain  the  (row)  rank  of  A,  pivoting  must  be  included  if  a  pivot 
element  is  0.  The  pivot  row  can  be  interchanged  with  a  row  of  the  active  matrix  with  a 
nonzero  entry  in  the  pivot  coiiimn.  If  no  such  exists  then  the  elimination  proceeds  to  the 
next  column.  The  rank  is  given  by  the  number  of  nonzero  elements  on  the  diagonal  of 
the  output  matrix. 

Determinant.  The  determinant  of  a  matrix  with  elements  in  72  can  be  defined  using 
permutations  and  using  that  definition  it  is  at  least  feasible  [5]  to  define  the  idea  of  an 
inverse  matrix  when  72.  is  a  commutative  ring  with  a  multiplicative  identity  element  1. 
(The  existence  of  1  is  needed  in  order  even  to  define  an  identity  matrix.  Commutativity 
is  necessary  for  the  identity  det  AB  =  (det  A)  (detB)  to  hold.) 

Howe\^er  these  conditions  are  not  sufficient  for  the  computation  of  the  determinant 
using  GE.  The  division-free  algorithm  does  not  readily  yield  det  A.  The  general  ring 
equivalent  of  (13)  clearly  demands  division  in  order  to  extract  the  deteiminant.  The 
fraction-free  Algorithm  14  however  is  well-defined  for  any  unique  factorization  domain. 
LTD,  since  ‘‘division'’  by  a  common  factor  is  now  available. 
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Algorithm  18  Computation  of  det  A  by  fraction-free  GE  in  a  unique  factorization  do- 
main 

Input  n  X  n  matrix  A  with  entries  in  a  UFD  Ti 
Compute 

for  1  =  1  to  n  ~  1 
factor  an 
for  j  =  2*  +  1  to  n 

for  fc  =  i  +  1  to  n  . 

djk  •—  O'iiO'jk 
ajj  0 

if  2  >  2  then  (remo\al  of  common  factors) 

for  j  =  7*  -f  1  to  n 
for  A:  =  i  +  1  to  71 
factor  ajk 

replace  ajk  by  product  of  all  its  factors  except  those  of 

Output  det  A  =  Unn 

Remark  24.  The  factorization  of  an ,  ajk  and  the  replacement  of  the  latter  with  a  reduced 
product  is  the  equivalent  of  division  in  the  integer  version,  Algorithm  14. 

Remark  25,  Such  an  algorithm  as  this  could  be  used  for  example  in  the  polynomial  ring 
R[x]. 

Remark  26.  Pivotmg  is  needed  in  the  event  that  any  zero-pivot  is  encountered.  This 
can  be  done  in  the  usual  waj*  by  interchanging  the  pivot  row  with  any  suitable  active  row 
—  and  negating  the  final  determinant  value. 

For  a  noncommutative  division  ring  —  that  is  a  ring  with  all  the  properties  of  a  field 
except  that  multiplication  is  not  commutative  —  the  original  Algorithm  1  can  be  modified 
to  provide  a  determinant. 

Algorithm  19  Computation  of  det  A  by  GE  in  a  division  ring  Tl 

Input  n  X  n  matrix  A  with  entries  from  TZ 
Compute 

for  2  =  1  to  n  —  1 
for  j  =  z  +  1  to  n 
for  A:  =  2  +  1  to  n 

0>jk  •“  {P'ik ! 

aji  I—  0 

Output  det  A  =  nr=i 

Remark  27.  Again,  of  course,  some  pivotwg  must  be  included  in  the  event  that  any 
pivot  is  zero. 
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5.3.  Solvability  for  systems  of  equations.  The  first  question  to  be  addressed  is 
“For  what  types  of  ring  can  we  solve  such  systems?’  In  [5]  we  see  that  the  notion  of  an 
inverse  matrix  is  definable  for  any  commutative  ring  TZ  with  1.  However,  in  such  a  general 
setting,  the  definition  of  a  matrix  being  nonsingular  can  be  confused. 

The  usual  (though  not  computationally  useful)  cofactor  definitions  can  be  used  to 
obtain  both  det  A  and  a  matrbc  B  satisfying 

AB  =  BA  =  {detA)/  (22) 

and  if  det  A  is  a  unit  in  Tl  then  (det  B  is  the  inverse  of  a  A.  (Re^ll  thaX  a  unit  u  in 
TZ  is  an  element  which  has  a  multiplicative  inverse  u  ^  such  that  uu  =  u  u  —  1.)  In 
such  circumstances,  the  system  Ax  =  b  has  a  solution. 

We  already  see  that  the  notion  of  nonsingularity  is  open  to  many  interpretations.  In  a 
more  general  ring  still,  we  can  define  and  compute  the  rank  of  a  matrix  using  Algorithm 
R1  for  example  which  allows  nonsingularity  of  A  to  be  defined  by  its  having  full  rank. 
This  definition  applies  even  when  detA  is  not  well-defined.  However,  if  detA  is  well- 
defined  and  is  nonzero,  then  A  has  full  rank  and  so  there  is  no  confusion  between  these 
two  definitions  of  nonsingular  when  they  both  make  sense.  From  (22)  it  is  apparent  that 
nonsingular  and  invertible  need  not  be  equivalent  since  det  A  could  be  nonzero  but  not  a 
unit. 

We  have  already  seen  that  for  detA  to  be  computed  by  GE,  we  require  that  TZ  is  b 
UFD.  In  such  circumstances  (and  assuming  A  is  nonsingular)  it  follows  that  a  system 
Ax  =  b  has  a  solution  in  the  field  of  fractions  Q,  This  solution  would  be  computed  by 
Algorithm  15  modified  to  the  “arithmetic”  of  7^  and  Q  with  the  operation  reduce  being 
interpreted  as  removal  of  any  common  factors  using  the  unique  factorization  property 
rather  than  the  gcd.  If  7?,  is  a  Euclidean  domain,  then  (a  suitably  modified  form  of)  the 
Euclidean  algorithm  for  obtaining  gi*eatest  common  divisors  can  be  used. 

Unique  factorization  domains.  In  order  to  obtain  the  solution  of  a  system  in  the 
field  of  fractions  Q  of  a  UFD  IZ  we  can  use  the  following  approach. 

Algorithm  20  Solution  of  a  system  by  GE  in  the  field  of  fractions  Q  of  a  UFD  1Z 

Input  n  X  n  matrix  A  and  an  n-vector  b  with  entries  in  a  UFD  TZ 

(Fraction-free  forward  elimination  as  in  Algorithm  UFDdet 
modified  for  right-hand  side) 

Compute 

for  i  —  1  to  71—1 
factor  an 
for  j  =  7  -h  1  to  71 

hj  :=  “  ajxbi 

for  A:  =  7-  +  1  to  7?. 

ajk  UjiQ'ik 

aji  :=  0 

if  7  >  2  then  (removal  of  common  factors) 

for  j  =  7  -h  1  to  71 
factor  bj 

replace  bj  by  product  of  all  its  factors  except  those  of  ai-i^^i 
for  /j  =  7  -h  1  to  7? 
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factor  ajk 

replace  ajk  by  product  of  all  its  factors  except  those  of 
(Back  substitution  as  in  Algorithm  15  modified  for  TZ,  Q) 

Initialize  rational  solution  m  :=  0;  n  :=  1  (solution  will  be  x,  =  ^i/^i) 

Rationalize  right-hand  side  p  :=  b;  q  :=  1  (right-hand  side  hi  =  VHQi) 

Compute 

for  j  =  n  downto  1 

reduce (factor  both  arguments,  remove  all  common  factors) 
for  z  =  1  to  i  -  1 

Pi  :=  Pi  ♦  nj  -  Uij  *  lUj  *  Qi 
Qi  :=  Qi  *  nj 

reduce  [pi,  Qi] 

Output  solution  [m,  n] 

Integral  domains.  Using  the  division-free  Algorithm  11  for  the  forward  elimination 
would  allow  the  UFD  solution  to  be  extended  to  the  case  where  7^  is  a  general  integral 
domain.  Again  a  solution  will  exist  in  the  field  of  fractions  Q  provided  the  matrix  is 
nonsingular.  In  this  situation,  just  as  with  divisionless  integer  arithmetic,  the  problems 
of  range  growth  become  potentially  severe. 

Algorithm  21  Division-free  solution  of  a  system  by  GE  in  a  the  field  of  fractions  Q  of 
an  integral  domain  TZ 

Input  77  X  n  integer  matrbc  A  and  right-hand  side  b 
Compute  (Division-free  forward  elimination) 
for  7  =  1  to  7?  —  1 

for  J  =  7  -i-  1  to  71 
bj  :==  Q>iibj 
for  =  7  -h  1  to  77 

aji  :=  0 

Initialize  rational  solution  m  :=  0;  n  :=  1  (solution  will  be  x,  =  mi/Ui) 

Rationalize  right-hand  side  p  :=  b;  q  :=  1  (right-hand  side  hi  =  pi/qi) 

Compute  (“Rational”  back  substitution) 
for  j  =  77  downto  1 

m.j  :=  pj;  nj  :=  Qj  *  Ujj 
for  7  =  1  to  j  —  1 
Pi  •=  Pi  *  j  “ 

Qi  :=  Qi  »  nj 
OutpiLt  solution  [m,  n] 

Remark  28.  The  potential  for  range  growth  is  showii  very  cJeaily  by  the  absence  here 
of  any  reduction  of  the  fractions  generated  in  the  back  substitution  phase. 
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Solution  in  Just  as  with  the  integers  themselves,  it  may  be  that  a  j>articular 
system  has  a  solution  in  the  original  ring  H  without  the  need  to  go  to  the  field  of  fractions 
and  rational  “arithmetic”.  In  the  case  where  7^  is  a  UFD  the  “integer”  solution,  if  it 
exists,  would  be  delivered  since  the  various  reductions  and  factorizations  would  result  in 
all  denominators  being  1.  However,  if  we  know  in  advance  that  a  particular  system  has 
a  solution  in  Tl,  then  at  least  the  back  substitution  phase  can  be  simplified  since  there  is 

then  no  need  to  introduce  the  fractions  at  all. 

The  simplifications  to  the  above  algorithms  are  then  obvious  and  we  do  not  detail 

them  here. 


6.  The  impact  of  parallelism 

The  recent  expansion  in  parallel  computing  has  had  profound  effects  on  many  are^  of 
computation  -  but  perhaps  nowhere  has  this  impact  been  so  marked  as  in  numerical  linear 
algebra.  There  are  many  different  parallel  architectures  including  many  special  purpose 
processors  for  tasks  such  as  signal  processing.  In  this  section,  we  restrict  attention  to  two 
basic  types  of  parallelism:  vector  processing  and  array  processors. 

The  former  are  characterized  by  the  highly  optimized  vector  pipeline  processors  of  the 
Cray  XMP  series.  These  machines  also  typically  featured  some  “true”  parallelism  in  that 
several  such  processors  are  available  to  be  used  simultaneously  working  on  different  data. 
The  latter  class  of  (massively  parallel)  array  processors  is  characterized  by  machines  such 
as  the  Connection  Machines  and  the  MasPar  MP*1  and  MP-2  series. 

Each  individual  architecture  has  its  own  particular  properties  but  there  are  some  gen¬ 
eral  properties  which  can  usefully  be  discussed  here.  For  much  more  detail  on  the  use  of 
vector  machines  for  solution  of  linear  systems  see  [18].  For  further  details  on  parallel  algo¬ 
rithms  and  array  processors,  [17]  provides  an  excellent  introduction  with  special  reference 
to  linear  algebra  problems.  The  discussion  of  basic  linear  algebra  operations  in  [7]  does 
not  pay  particular  attention  to  any  one  architecture  but  does  address  architecture-related 
issues  through  discussions  of  saxpy  and  other  special  operations  and  their  relations  to 
storage  schemes.  An  introduction  to  some  of  these  ideas  and  to  performing  basic  linear 
algebra  operations  on  an  array  processor  can  also  be  found  in  [3]. 

0,1,  Vector  and  pipeline  architectures.  By  a  “vector”  architecture,  we  mean  a 
computer  which  is  designed  to  perform  operations  on  vector  quantities  at  high  speed.  This 
is  most  frequently  achieved  through  efficient  use  of  a  pipeline  processor.  In  such  a  machine 
the  operation  of  adding  tw^o  vectors  for  example  is  greatly  enhanced  by  allowing  certain 
subtasks  (such  as  reading  form  memory,  adding,  rounding  and  normalizing,  and  writing 
to  memory)  to  be  performed  simultaneously.  From  the  point  of  view*  of  programming  a 
vector  machine,  it  is  usually  valid  to  think  of  the  operation 
for  ?‘  =  1  to  n 
C,  Uj  “4“ 

as  a  single  vector  operation 

c  :=  a-hb 

w*here  each  component  of  c  is  computed  at  the  same  time. 

There  are  some  occasions  when  it  is  necessary  to  pay  more  careful  attention  to  the 
pipelined  nature  of  such  a  vector  operation  but  for  the  most  part  thinking  of  it  as  a  vector 
operation  is  appropriate.  In  fact  most  ve<?tor  machines  are  designed  to  peifoim  a  typical 
saxpy  operation  such  as  v  :=  ax  -f  y  (“scalar  a  x  plus  y  )  very  efficiently. 
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Note  that  the  innermost  loop  of  the  basic  ijk  form  of  Gauss  elimination,  Algorithm  1, 
is  performing  a  row-saxpy.  this  implementation  of  GE  would  therefore  be  well-suited  to  a 
vector  computer  in  which  matrices  are  stored  by  rows.  Similarly  the  ikj  form,  Algorithm 
5,  also  has  a  saxpy  at  its  heart  -  but  this  time  a  column-saxpy  making  this  version  of  GE 
better  suited  to  a  vector  processor  with  matrices  stored  by  columns. 

To  see  the  real  effect  of  the  vector  processing  here,  we  consider  the  column-oriented 
version.  For  notational  convenience  we  denote  that  part  of  the  i-th  column  of  A  below  the 
rpnin  diagonal  by  c(9  so  that  the  fc-th  component  of  this  column  vector  is  c]^’  =  ai+k,i- 
The  part  of  the  right-hand  side  vector  below  the  i-th  position  is  denoted  b^’).  With  this 
notation  the  vector  form  of  Algorithm  5  is 

Algorithm  22  Basic  GE:  Column  vector  form 

Input  nxn  matrix  A  (and  right-hand  side  b  if  solving  a  sj'stem) 

Compute 

for  i  =  1  to  n  —  1 
mb)  ;=  cb)/aij 

ab)  :=  0 

bb)  :=  bb)  -  b^mb)  (if  solving  a  system) 
for  k  =  i  +  1  to  n 

:=  —  au-mb) 

Output  (modified)  matrfac  A  (and  b  if  solving  a  system) 

The  major  impact  on  the  complexity  of  the  algorithm  is  immediately  apparent:  the 
replacement  of  the  y-loop  by  a  single  vector  operation  reduces  the  overall  complexity 
to  O(n^).  Computation  of  the  multiplier  vector  is  consists  of  division  of  a  vector  by  a 
constant  which  can  be  achieved  with  a  single  (scalar)  reciprocation  followed  by  a  saxpy 
in  which  one  of  the  vectors  is  0.  Similarly  the  update  of  the  right-hand  side  vector  is  a 
column-saxpy. 

The  vector  operation  counts  for  this  algorithm  are  summarized  in  Table  6  along  with 
those  for  the  bade  substitution  phase  which  is  a  similarly  vectorized  version  of  Algorithm 
8- 

TABLE  6  Vector  floating-point  operation  counts  for  column  GE  solution  of  an  n  x  n 
linear  system  Ax  =  b. 


Operation 

Forward  elimination 

Back  substitution  TOTAL 

Scalar  / 
saxpy 

n  —  1 

1  („!)(„ +  2) 

n  (2n  - 1) 

n-1  L(„  _  i)(n-(-4) 

In  the  case  of  a  pipeline  computer  as  opposed  to  a  truly  parallel  vector  machine,  some 
of  the  vector  lengths  would  become  too  short  to  justify  using  the  pipeline.  Nonetheless  we 
see  that  the  potential  saving  achieved  by  even  the  most  elementary  parallelism  is  great. 
In  the  case  of  a  “true’’  vector  processor,  we  have  made  an  implicit  assumption  that  the 
limit  on  the  vector  length  is  at  least  7i. 

We  see  that  adding  a  “vector  dimension”  to  the  processor  power  has  reduced  the 
algorithmic  complexity  by  an  order  of  magnitude.  In  the  case  of  an  ar  ray  processor  with 
at  least  processors  we  might  therefore  hope  to  reduce  the  Gauss  elimination  algorithm 
to  O  (?t)  parallel  floating-point  operations,  or  paraftops. 
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6  2  Array  architectures.  In  this  section,  we  consider  the  use  of  massively  parallel 
a^rav  processors  for  Gauss  elimination.  Such  a  computer  architecture  consi^s  of  a  (typi¬ 
cally  square  or  rectangular)  array  of  procesing  elements  operating  in  the  SIMD  paradigm. 
That  is  to  sav  that,  at  any  given  time,  each  processor  is  either  performing  the  same  in¬ 
struction  on  its  own  data  or  is  idle.  (SIMD  is  an  acronym  for  Single  Instruction  Multiple 

ThT^vidual  processors  in  such  an  array  are  usually  much  less  powerful  than  would 
be  encountered  on  a  serial  machine.  The  power  of  the  paralleUsm  comes  from  its  scale. 
Even  if  individual  operations  are  much  slower,  adding  two  matrices  using  a  single  paratlop 
with  the  instruction 

C  A B 

is  likely  to  be  significantly  faster  than  the  conventional  double  loop 
for  7  =  1  to  71 
for  j  =  1  to  71 

dj  :=  dij  +  bij 

with  its  71^  flops  —  especially  when  n  gets  large.  r  ^  i  4.  Ar\c\r 

T>'pical  parallel  arrays  have  at  least  1024  =  32  x  32  processors.  Arrays  of  at  least  4096 
are  now  more  common.  On  such  an  array  the  comparison  for  matrix  addition  k  therefore 
between  one  (slow)  paraflop  and  4096  (fast)  serial  floating-point  operations.  Even  if  the 
individual  arithmetic  operations  are  slower  by  a  factor  of  32,  the  overall  performance 
would  be  speeded  up  by  a  factor  of  128  which  is  certainly  worth  achieving.  Many  linear 
algebra  operations  can  be  readily  modified  to  be  performed  on  an  array  processor  with 
close  to  the  maximal  speed-up.  Matrix  multiplication,  for  example  can  be  reduced  to  just 
Oln)  paraflops  on  an  ti  x  ti  array. 

For  simplicity,  we  shall  suppose  throughout  this  section  that  the  matrices  under  con¬ 
sideration  fit  the  array  —  so  that  both  are  n  x  n  arrays.  For  larger  matrices,  the  sort  of 
algorithms  discussed  here  can  be  modified  to  deal  with  the  matrices  in  pieces.  For  smaller 
matrices,  the  algorithms  described  can  be  used  but  will  not  necessarily  take  full  advantage 

of  the  processor  array.  . 

The  biggest  new  cost  associated  with  implementing  algorithms  such  as  Gauss  elimi¬ 
nation  or  LU  factorization  on  an  array  processor  comes  form  the  inter-processor  commu¬ 
nications  which  are  needed.  This  cost  must  be  weighed  alongside  all  the  arithmetic  and 
other  costs  in  deciding  what  architecture  should  be  used  for  a  particular  problem.  In  this 
section  we  are  concerned  only  with  implementing  GE  on  such  a  processor. 

In  keeping  with  the  notation  of  the  previous  section,  we  shall  denote  by  r  ‘  the  !-th  row 
of  the  current  matrix  .4.  We  shall  also  define  two  logical  arrays  rows  and  columns  which 
are  used  to  control  which  processors  perform  particular  operations.  For  the  t-th  stage  of 
GE,  these  sets  will  are  rows  =  {j  :  j  >  ?}  and  columns  =  {k:k>i}.  The  “active’’  part 
of  the  matrix  would  then  be  identified  with  the  mask 

if  (rows)  and  (columns)  ...  ,  i  ^  ■  ij 

which  has  the  effect  that  only  those  processors  for  which  both  j  >  t  and  k  >  i  would 

perform  the  instructions  covered  by  this  test.  .  .  , 

In  Algorithm  23.  below,  for  the  case  of  linear  systems  we  consider  the  situation  of 

solving  a  matrix  system 

.4A'  =  B 

where  B  can  have  as  many  as  n  columns.  There  is  no  additional  cost  associated  with  this. 
We  shall  denote  the  r-th  row  of  the  current  matrix  B  by  b^'^  (Note  this  is  different  from 


NAWCADPAX-96- 1 94-TR 


Gauss  Elimination:  Workhorse  of  Linear  Algebra  39 


the  notation  for  Algorithm  22  for  a  vector  processor  which  was  column-oriented.)  The 
notation  used  in  the  last  section  is  retained. 

The  other  point  to  be  emphasized  before  detailing  the  algorithm  is  that  arithmetic 
operations  referred  to  here  are  array  operations  not  matrix  operations.  Thus  the  product 
MA  refers  to  the  elementwise  product  of  these  two  arrays. 

Algorithm  23  Basic  GE:  Array  processor  form 

Input  nxn  matrix  A  (and  right-hand  side  B  if  solving  a  system) 
Compute 

for  i  =  1  to  n  —  1 

rows  :=  {j  :  j  >  r} 
columns  :=  {k:  k>  i} 

copy  to  all  rows  to  generate  the  array  rowa 

copy  to  all  rows  to  generate  the  array  rowb 

copy  to  all  columns  to  generate  the  array  M 
if  (roic’s) 

B  :=  jB  -  M  *  rowb  (if  solving  a  system) 
if  {columns)  A::=^  A--  M  *  rowa 
Output  (modified)  matrix  A  (and  B  if  solving  a  system) 

Remark  29.  The  arithmetic  complexity  of  this  algorithm  is  easy  to  compute:  There  are 
just  5  paraBops  used  at  each  stage  so  that  the  total  count  is  5  (n  —  1)  paraflops  consisting*  of 
(n  -  1)  parallel  divisions,  2  (n  -  1)  parallel  multiplications  and  2  (n  -  1)  parallel  additions. 

It  is  also  easy  to  see  that  Algorithm  23  is  wasteful.  At  stage  i,  all  processors  in 
the  first  i  rows  are  inactive.  However,  at  no  additional  cost,  the  elimination  could  be 
performed  above  the  diagonal  simultaneous  with  the  triangular  factorization.  That  is  we 
could  perform  the  Gauss- Jordan  algorithm  for  the  same  cost  as  GE.  (Recall  that  for  a 
serial  machine,  Gauss- Jordan  is  twice  as  expensive  eis  GE.)  The  effect  of  this  is  that  the 
back  substitution  is  replaced  by  solving  a  diagonal  system  -  a  very  simple  task  -  especially 
on  a  parallel  computer.  The  only  change  needed  is  to  modify  the  definition  of  the  row 
mask  to  rows  :=  {j  :  j  ^  i}- 

The  solution  of  the  resulting  diagonal  system  would  be  achieved  by  simply  copying  the 
diagonal  entries  across  their  respective  rows  and  then  using  the  single  parallel  division: 

A'  :=  Bfdiag 

This  even  raises  the  possibility  that  there  may  be  circumstances  in  which  computing 
a  matrix  inverse  may  be  a  practical  and  desirable  computation.  Of  course,  it  should  be 
recalled  that  the  Gauss-Jordan  algorithm  is  numerically  less  stable  than  Gauss  elimina¬ 
tion.  For  ill-conditioned  matrices,  therefore  it  may  be  better  not  to  use  this  variation  - 
however  for  such  matrices  we  would  probably  want  to  use  a  more  stable  algorithm  than 
GE,  too! 

Pivoting.  As  in  all  other  cases  of  GE,  it  is  possible  that  zero  or  small  pivots  will  be 
encountered  and  that  pivoting  is  needed.  In  most  parallel  programming  languages  there 
is  a  built-in  reduction  algorithm  for  locating  the  maximum  element  in  an  array.  Therefore 
finding  the  conventional  pivot  element  for  GE  is  a  simple  task.  The  choice  between  row 
interchanges  and  storing  a  pivot  vector  is  much  the  same  as  for  serial  implementations. 
The  pivot  vector  could  easily  be  replaced  in  this  context  with  another  mask  us€d  so  that 
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the  search  for  the  pivot  element  would  only  consider  rows  which  had  not  yet  been  used 
as  a  pivot  row. 

6.3.  Parallelism  for  long  integer  arithmetic.  As  a  final  comment  on  the  impact  of 
parallelism  on  performing  Gauss  elimination,  we  return  to  the  question  of  integer  compu¬ 
tation  and  the  need  for  long  integer  arithmetic  to  accommodate  the  growT.h  in  the  dynamic 
range.  In  Section  3.7  we  considered  the  requirements  of  such  long  integer  arithmetic  and, 
in  particular,  the  convolution  form  of  the  product  of  such  integers. 

The  problem  was  described  as  finding  the  product  of 

K-\  A-i 

m  =  ^  diR'  ■  ^ 

»=o 

where  R  is  the  effective  radix  resulting  form  breaking  a  long  wordlength  integer  into 
shorter  integer  pieces.  Recall  that  in  (17),  the  product  is  given  by 


m*n  — 


which  is  a  convolution  product  of  the  coefficient  vectors.  Each  coefficient  product  aiPk-i 
is  a  base-R*  digit  which  we  wite  as  =  lk,iR  +  where  each  7fc,i  and  is  a 

base-R  digit.  The  product  (17)  can  therefore  be  written  as 


m  * 


n  =  +  4,i  -  13 


+  6k,i  R' 


where  ')k,i  —  Sk,i  =  0  for  f  >  fc. 

An  array  processor  could  be  used  to  accelerate  this  algorithm  greatly.  Firstly,  all  the 
coefficient  products  in  the  inner  sum  in  (17)  can  be  computed  simultaneously.  The  re¬ 
alignment  of  the  upper  and  lower  halves  of  these  products  needed  for  the  inner  summation 
in  the  last  equation  is  then  a  simple  shift  of  data  to  a  neighboring  processor  and  these 
sums  could  then  also  be  performed  simultaneously.  The  final  carries  would  be  the  only 

part  that  requires  serial  processing.  ,  •  • 

These  arithmetic-algorithmic  considerations  are  analogous  to  the  design  decisions 
which  are  used  with  R  =  2  in  the  design  of  hardware  arithmetic  units  where  bit-parallelism 
is  exploited  wherever  possible. 


7.  Noisy  matrices 

In  this  section  we  address  the  problem  of  computing  with  “noisy”  matrices.  That  is  we 
shall  assume  that  the  underlying  matrbc  A  is  contaminated  with  a  noise  matrbc  E  so  that 
the  matrix  with  which  we  must  compute  is 

A  =  A  +  E  (23) 

We  shall  also  assume  here  that  the  level  of  the  noise  can  be  estimated  so  that  we  have 
bounds  on  the  elements  of  E  such  as 


\eij\  <  S 


(24) 
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or  that  we  know  the  elements  of  E  are  normally  distributed  random  variables  with  mean 
0  and  standard  deviation  o.  In  the  latter  case,  it  follows  that  for  each  i*j 

\eij\<2a  (25) 

with  probability  about  95%.  (See  any  standard  text  on  Statistics  such  as  [8].)  FVom 
the  similarity  in  (24)  and  (25)  it  appears  that  the  Gaussian  noise  situation  can  be  well- 
modeled  by  the  uniform  noise  model  at  least  in  terms  of  trying  to  estimate  error  bounds 
or  tolerances.  For  the  remainder  of  this  section,  we  use  the  uniform  noise  model  (24). 

Provided  that  S  is  not  too  large,  the  effect  of  using  A  rather  than  A  is  covered  by 
the  standard  error  analyses  for  all  the  basic  linear  algebra  problems  with  matrices  of  full 
rank.  Essentially  S  becomes  a  bound  on  the  absolute  errors  of  the  original  data  or  E  plays 
the  role  of  6 A  in  the  error  analysis  described  in  Section  2.5.  (See  [3]  or,  for  an  extensive 
treatment,  [7]  for  more  details  of  this  error  analysis.)  One  situation  in  which  the  noisy 
matrbc  (23)  can  be  expected  to  yield  substantially  different  results  from  A  is  when  A  is 
rank  deficient.  The  determination  of  the  “effective”  rank  of  rank  deficient  matrices  is  an 
important  aspect  of  several  naval  applications. 

7.1.  Determination  of  effective  rank.  By  deteipination  of  the  effective  rank,  we 
mean  finding  rank(>l)  even  though  we  are  given  only  A  and  the  bound  (24)  on  the  noise 
matrbc.  We  assume  here  that  all  matrices  have  dimension  nxn  and  that  A  has  rank  less 
than  n.  Typically,  we  would  expect  A  to  be  of  full  rank  since  the  noise  matrbc  is  unlikely 
to  reflect  the  same  linear  dependencies  among  its  rows  (or  columns)  as  those  of  A. 

In  the  no-noise  situation,  Algorithm  10  delivers  the  rank  of  a  matrbc  simply  by  counting 
the  number  of  nonzero  elements  on  the  diagonal  of  the  upper  triangular  factor  U  result¬ 
ing  from  LU  factorization  with  partial  pivoting.  Of  course,  even  roundoff  errors  mean 
that  testing  these  diagonal  elements  against  zero  is  inadequate.  For  very  ill-conditioned 
matrices,  this  effect  could  be  severe.  But  for  other  matrices,  we  may  expect  that  these 
zeros  will  be  replaced  with  small  values  and  that  the  true  nonzero  diagonal  entries  are 
substantially  larger. 

In  [12],  it  is  stated  that  GE  (or  LU  factorization)  cannot  satisfactorily  distinguish 
between  the  “true”  nonzeros  and  the  “should-be- zeros” .  However  this  conclusion  a]> 
pears  to  be  based  on  an  example  in  which  a  divisionless  version  of  Gauss  elimination  is 
applied  to  a  small  example.  If  Algorithm  10  (even  without  pivoting)  is  applied  to  this 
example  the  results  are  very  clear.  In  their  application  the  diagonal  entries  are  given 
as  2.998,  -12.922,0.577  which  should  be  compared  with  singular  values  10.331,  2.644  and 
0.0064.  However,  with  or  without  pivoting,  conventional  GE  yields  a  diagonal  with  entries 
2.998,4.315,-0.014  which  shows  a  much  clearer  break  between  the  “large”  values  and  the 
small  one.  In  this  particular  example,  S  =  0.01  so  that  the  “zero”  has  been  preserved  to 
the  same  order  of  magnitude  as  S. 

The  question  is  how  well  can  GE  do  as  a  rank  determination  algorithm  in  the  presence 
of  noise?  Clearly  one  3x3  example  does  not  justify  its  general  applicability.  Some 
experiments  were  conducted  using  MATLAB  with  randomly  generated  noisy  matrices 
where  the  bound  5  was  small  relative  to  the  range  of  values  for  the  elements  of  A. 

Generation  of  random  noisy  matrices.  In  order  to  generate  random  noisy  matri¬ 
ces  with  known  effective  rank  t\  first  a  set  of  r  linearly  independent  vectors  with  random 
entries  in  a  specified  interval  v\7is  generated.  Next  random  linear  combinations  of  these 
row  vectors  were  used  to  genei'ate  the  rows  of  the  matrix  A .  The  coefficients  of  these  linear 
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combinations  were  also  from  a  specified  range.  Finally  a  complete  n  x  n  noise  matrix  with 
random  elements  in  the  interval  [-5,5]  was  added  to  yield  A. 

Of  course  there  is  a  remote  possibility  that  the  underlying  matrix  A  could  have  rank 
less  than  r.  This  never  happened.  An  example  of  the  MATLAB  code  used  for  this  is. 

function  A  =  testmat(n.r): 

a=20*rand(n)-10: 

A=zeros(n): 

for  i=l:n 
for  j=l:r 

A(i,:)=A(i,:)-f(6*rand-3)*a(j,-.): 

end 

end 

function  E=noisemat(n,s): 

E=s*(2*rand(n)-1); 

Ahat=testmat(  n ,  r)-rnoisemat(  n  ,s) ; 

The  tolersince  level.  With  partial  pivoting,  all  multipliers  used  in  the  LU  factor¬ 
ization  are  necessarily  less  than  1.  Assuming  that  the  noise  is  significantly  greater  than 
the  roundoff  error  effect  but  significantly  smaller  than  the  range  of  the  matrix  entries  ,  it 
follows  that  its  effect  on  the  LU  factors  is  bounded  by  nS  and  so  this  could  be  used  as  a 
likely  tolerance  for  determining  the  effective  rank  of  the  matrfac. 

The  MATLAB  experiments  that  were  performed  suggest  that  this  is  indeed  a  satis¬ 
factory  tolerance  to  use  most  of  the  time.  Some  experiments  were  carried  out  using  the 
MATLAB  LU  factorization  routine.  The  initial  tests  were  based  on  perturbations  of  a 
7x7  matrix  which  was  used  as  the  basis  of  the  analysis  in  [12].  This  matrix  is 

■  3  2  1  7  4  5  3  ■ 

1  4  2  6  5  10  3 

8  1  5  13  5  7  0 

A=  4  2  7  15  11  11  4  (26) 

1  2  1  3  2  5  1 

2  1  3  5  3  5  0 

3  10  1  5  2  21  1  _ 

which  has  true  rank  4. 

The  MATLAB  function  randn  was  used  to  generate  noise  matrices  whose  elements 
were  normallv  distributed  and  with  specified  standard  deviation.  The  first  phase  of  the 
experiment  was  to  obtain  the  effective  rank  of  many  such  perturbations  by  counting  the 
number  of  diagonal  entries  in  U  which  exceeded  some  threshold  absolute  \^lue.  Two 
hundred  perturbations  were  used  in  each  test.  The  value  of  the  standard  deviation  and 

the  tolerance  used  were  varied  between  runs. 

It  quickly  became  apparent  that  choosing  a  tolerance  close  to  ncr  provided  the  b«t 
results.  The  following  table  of  results  for  a  =  0.3  illustrate  this  point. 
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TABLE  7  Rank  tests  using  LU  factorization.  7x7  matrix  (26).  a  —  0.3 

Shows  number  of  times  each  rank  was  obtained  with  200  perturbations  for  each  tolerance 


Tolerance 

Rank=3 

Rank=4 

Rank=5 

Rank=6 

Rank— 7 

a 

0 

0 

10 

72 

118 

1 

0 

60 

97 

39 

4  . 

ba 

0 

155 

43 

2 

0 

lOc 

50 

150 

0 

0 

0 

When  the  test  was  repeated  for  <7  =  0.1, 1000  perturbations  using  tolerance  lOa  yielded 
rank  4  996  times  and  rank  5  just  4  times.  With  a  smaller  x'alue  of  o  the  results  were  ®ven 
more  successful.  For  both  these  small  noise  levels,  using  to\  =  ha  gave  approximately  75% 
success.  These  exi>eriments  suggest  that  GE  with  the  appropriate  tolerance  could  be  the 
basis  of  a  successful  rank  determination  algorithm  when  the  noise  level  is  known.  It  is 
also  already  apparent  that  finding  the  right  tolerance  to  use  is  of  critical  importance  and 
the  results  may  be  highly  sensitive  to  this  choice. 

One  way  of  reducing  this  effect  is  to  use  the  LU  factorization  as  a  basis  for  a  proba¬ 
bilistic  approach  to  the  problem.  To  this  end  a  second  set  of  tests  was  run  in  which  the 
original  matrix  A  is  perturbed.  The  LU  rank  test  used  in  the  previous  experiments  is 
then  applied  to  a  number  of  perturbations  of  this  matrix.  The  effective  rank  would  then 
be  determined  by  the  most  frequent  “rank"  for  these  and  some  estimate  of  the  confidence 
we  should  have  in  the  result  is  provided  by  the  complete  array  of  results.  Thus  a  single 
test  would  produce  results  similar  to  a  row  of  Table  7. 

In  the  experiments  200  perturbations  w-ere  again  used.  Clearly  this  would  be  too 
exi>ensive  since  for  even  a  poorly  conditioned  example  we  would  expect  to  be  able  to 
complete  the  singular  value  decomposition  of  A  with  much  less  effort  than  200  LU  factor¬ 
izations.  The  reeison  for  using  so  many  test  runs  here  is  to  try  to  eliminate  falsely  positive 
results  based  on  small  samples. 

The  tests  were  performed  using  three  different  noise  levels  a  —  10"  ,0.1, 0.3  «»ch  wdth 
two  tolerance  levels  7a  and  lOo.  The  combined  results  of  four  runs  of  each  case  are 
summarized  in  Table  8.  Each  run  uses  a  different  underlying  perturbation  of  the  basic 
7x7  matrix. 

TABLE  8  R.ank  tests  using  LU  factorization.  Perturbations  of  the  7x7  matrix  (26). 

Table  shows  percentage  for  each  rank  was  obtained  in  each  case  and  the  resulting 

estimated  rank. 


Std  Deviation 

Tolerance 

Rank=3 

Rank=4 

Rank=5 

Rank=6 

Est.  Rank 

7a 

0 

71.0 

28.0 

1.0 

4 

0.1 

la 

0 

76.25 

23.38 

0.38 

4 

0.3 

la 

2.88 

75.25 

20.5 

1.38 

4 

10- ® 

10a 

0 

91.5 

8.5 

0 

4 

0.1 

10a 

0 

92.63 

7.-38 

0 

4 

0.3 

10a 

5.5 

88.13 

6.38 

0 

4 
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This  evidence  suggests  that  for  a  matrix  of  this  dimension  a  tolerance  of  lOcr  is  likely 
to  prove  very  reliable  in  predicting  the  effective  rank  of  a  matrix  from  the  LU  factorization 
of  even  a  small  number  of  perturbations  of  the  original  matrix.  As  would  be  expected  the 
performance  appears  to  be  beginning  to  deteriorate  as  the  noise  level  rises. 

The  biggest  drawback  to  using  GE  in  this  context  is  that  setting  the  appropriate 
tolerance  is  difficult  and  especially  as  the  dimension  of  the  matrices  increases  the  cost  of 
performing  several  LU  factorizations  of  perturbed  copies  may  be  too  expensive.  Other 
approaches  to  this  problem  may  be  more  successful  for  higher  dimensions. 

Limited  further  experiments  were  conducted  with  matrices  of  higher  dimension  and 
without  the  probabilistic  element.  Specifically  within  the  context  of  low  rank  matrices  of 
larger  dimension,  the  instability  of  LU  factorization  becomes  apparent.  On  experiments 
using  randomly  generated  20  x  20  matrices  with  rank  4,  there  were  examples  in  which 
the  effective  rank  was  correctly  predicted  by  the  LU  factorization.  These  typically  had 

moderate  condition  numbers.  n  _ 

One  particular  example  generated  using  the  functions  described  above  with  n  =  20,  r  - 
4,  s  =  0.5  had  condition  number  reported  by  MATLAB  as  2.4  x  10®.  The  LU  factorization 
using  a  tolerance  of  ns  yielded  an  effective  rank  of  just  2.  (The  SVTD  gave  the  correct 
result  as  did  a  least  squares  based  algorithm.)  Studying  the  diagonal  entries  of  U  it  was 
apparent  that  there  was  no  readily  recognizable  alternative  tolerance  w^  available.  The 
two  “large"  entries  were  approximately  -39  and  31,  the  next  four  were  7.0,  6.8,  5.7  and 
4.5  followed  by  three  between  2  and  2.5  and  all  except  one  were  greater  than  0.5. 

Trying  repeated  perturbations  of  the  matrix  as  in  the  second  phase  of  the  above 
experiments  did  not  provide  any  clear  conclusion  either.  After  28  such  perturbations  had 
been  tested,  the  results  indicated  rank=4  17  times,  rank=3  10  tim^  and  rank=2  once. 
With  a  mean  rank  estimate  of  3.57  this  is  clearly  not  conclusive.  Note  too  that  28  LU 
factorizations  of  different  20  x  20  matrices  required  approximately  140,000  flops  which 
does  not  compare  favorably  with  the  31,000  required  for  the  S\T)  for  this  same  matrix. 

Gauss  elimination  may  be  the  answer  to  many  linear  algebra  problems  but  effective 
rank  determination  for  low  rank  matrices  of  even  moderately  large  dimension  in  the  pres¬ 
ence  of  noise  is  probably  best  approached  by  alternative  techniques. 

7.2.  Other  approaches.  The  most  reliable  and  commonly  used  technique  for  this 
problem  is  based  on  the  SVD.  This  algorithm  too  is  sensitive  to  the  tolerance  level  and 
that  is  the  main  point  of  [12].  However  in  general  the  singular  value  decomposition  is 
usually  expensive  to  compute. 

Tlie  work  of  [6]  is  based  on  an  alternative  which  makes  use  of  the  coefficients  of  the 
characteristic  polynomial  of  the  matrbc.  This  is  a  statistical  approach  which  relies  on  a 
finite  algebraic  algorithm  -  but  one  which  is  also  likely  to  be  expensi\e  computational!)’ 
since  it  is  most  naturally  suited  to  working  from  full  rank  downwards  until  the  effective 
rank  is  revealed.  The  problem  here  is  that  the  full  rank  determination  requires  the  com¬ 
putation  of  the  determinant  and  so  is  already  as  expensive  as  the  LU  factorization.  There 
may  be  alternative  arrangements  of  the  algorithm  which  would  overcome  this  difficult). 

A  further  alternative  currently  being  investigated  uses  least  squares  appi-oximation 
[24].  The  idea  here  is  to  have  an  algorithm  which  works  upwards  from  rank  1  to  determine 
the  rank  of  a  matrix  which  is  known  to  be  of  low  rank  more  quickly.  The  basic  idea  is 
to  approximate  each  row  of  the  matrix,  in  a  least  squares  sense,  witli  a  multiple  of  the 
largest  row.  If  the  residual  matrix  is  small  enough  then  the  effective  rank  is  1.  Otherwise 
the  process  is  repeated  with  the  remaining  rows.  The  early  results  with  this  approach  are 


NAWCADPAX-96-194-TR 


Gauss  Elimination:  Workhorse  of  Linear  Algebra 


45 


encouraging  and  the  investigation  is  continuing.  For  the  specific  example  quoted  above 

the  correct  rank  was  determined  using  just  8000  flops  representing  a  saving  of  nearly  75% 

compared  with  the  S\T>based  algorithm. 
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